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SECTION  1 


INTRODUCTION 

1-1  SCOPE  AND  BACKGROUND. 

This  report  documents  the  results  of  an  investigation 
conducted  with  the  objectives  of  (1)  developing  an  improved 
analytical  model  for  calculating  closure  and  failure  resistance 
of  deep-based  tunnels  and  (2)  relating  the  improved  models  to 
results  of  laboratory  studies. 

This  study  is  a  continuation  of  "An  Investigation  of 
the  Failure  Resistance  of  Rockbolted  Tunnels  for  Deep  Missile 
Basing"  (AA,  1983).  That  report  described  an  improved  design 
methodology  for  minimally  hardened  tunnels  that  relies  on 
elastoplastic  models  of  a  cylindrical  tunnel  in  an  infinite 
homogeneous  rock  mass.  In  these  models,  the  strength  of  rock  is 
simply  defined  in  terms  of  a  friction  angle  and  a  cohesion,  and 
the  volumetric  expansion  of  rock  during  failure  in  terms  of  a 
dilatancy  angle  (the  dilatancy  angle  controls  the  rate  of 
increase  of  the  inelastic  volume  change  of  the  material  with  the 
plastic  shear  strain).  With  respect  to  previous  methodologies 
(see  Reed  et  al.,  1983,  for  a  state-of-the-art  review),  three 
significant  improvements  were  achieved  in  the  course  of  that 
investigation : 

•  Consideration  of  an  arbitrary  dilatancy  angle 
with  values  anywhere  between  zero  (no  inelastic 
dilatation)  and  the  value  of  the  friction  angle 
(the  maximiun  theoretical  limit),  thus  providing  a 
generalization  of  the  models  of  Newmark  (1970) 
and  Hendron  and  Aiyer  (1971)  for  hydrostatic 
far-field  stress  loading. 

•  Consideration  of  a  deviatoric  component  in  the 
far-field  stress  (previous  analytical  models  were 
restricted  to  consideration  of  a  hydrostatic 
far-field  loading),  with  the  consequence  that  the 
model  can  predict  ovalling  of  the  tunnel  during 
closure. 


•  Development  of  design  charts  based  on  the  semi- 
analytical  nonhydrostatic  model.  These  charts 
constitute  a  powerful  and  quick  means  of  calcu¬ 
lating  the  support  pressure  to  be  provided  by  a 
support  system,  to  limit  closure  to  within  a 
preset  amount. 

That  investigation  concluded  that  prediction  of  tunnel 
closure  (and  thus  of  the  support  pressure)  is  very  sensitive  to 
the  assumed  (constant)  value  of  the  dilatancy  angle,  and  that 
the  dilatation  model  of  rock  based  on  a  constant  value  of  the 
dilatancy  angle  is  inadequate  for  this  class  of  problems.  In 
practice,  the  rate  of  increase  of  the  plastic  volume  change  with 
the  plastic  shear  strain  has  not  been  observed  to  be  constant  in 
rocks;  it  changes  with  the  aunount  of  plastic  deformation 
approaching  zero  as  the  "damage"  increases.  In  contrast,  the 
theoretical  models  based  on  the  assumption  of  a  constant  dila¬ 
tancy  angle  predict  that  there  is  no  upper  limit  on  the  maximum 
inelastic  volume  increase  that  the  material  can  experience.  The 
inadequacy  of  the  assumption  of  a  constant  dilatancy  angle  is 
particularly  severe  in  the  tunnel  problem,  which  is  charac¬ 
terized  by  a  high  distortional  strain  field  in  the  rock  mass. 

The  principal  objective  of  the  present  investigation 
was  then  to  develop  an  elastoplastic  model  of  a  tunnel  under 
either  hydrostatic  or  nonhydrostatic  loading,  using  a  more 
realistic  dilatation  model  for  the  rock.  The  improved  dilata¬ 
tion  model  considered  in  this  report  assumes  an  exponential 
decay  of  the  dilatancy  angle,  from  an  initial  value  equal  to  the 
friction  angle.  It  is  a  simple  law  that  has  the  advantage  of 
depending  on  only  one  physically  meaningful  parameter;  the 
maximum  plastic  volume  change  that  the  material  can  experience. 
Secondary  objectives  of  this  investigation  consisted  of  a  review 
of  the  rock  dilatation  phenomenon  and  an  analysis  of  physical 
model  experiments,  with  the  purpose  of  validating  the  improved 
dilatation  model. 


1-2  REPORT  ORGANIZATION. 

The  main  body  of  this  report  is  comprised  of  three 
major  sections  covering  respectively  (1)  the  phenomenon  of  rock 
dilatancy,  (2)  a  general  description  of  the  elastoplastic  models 
of  a  tunnel  under  hydrostatic  and  nonhydrostatic  loading,  and 
(3)  an  analysis  of  physical  model  tests.  Mathematical  deriva¬ 
tion  and  description  of  numerical  procedures  have  been  docu¬ 
mented  in  three  Appendices. 
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SECTION  2 


DILATANCY  OF  ROCKS 

2-1  CAUSES  OF  DILATANCY. 

The  phenomenon  of  volumetric  expansion,  or  dilatancy, 
has  long  been  observed  during  shear  deformation  of  densely 
packed  granular  media,  where  it  is  associated  with  relative 
movement  of  grains  and  is  a  geometrical  necessity  for  deforma¬ 
tion  to  occur.  Dilatant  behavior  of  rocks  during  failure  was 
observed  first  by  Bridgman  (1949)  during  compression  tests  on 
soapstone  and  calcite  marble.  Dilatancy  during  uniaxial  and 
triaxial  compression  tests  was  subsequently  confirmed  for  a 
large  variety  of  rocks:  norite  and  quartzite  (Bieniawski, 
1967),  granite  (Brace  et  al.,  1966;  Zoback  and  Byerlee,  1975), 
marble,  sandstone,  limestone,  etc.  Dilatancy  therefore  appears 
to  be  a  pervasive  property  of  many  rocks. 

A  typical  result  of  a  triaxial  compression  test  on  a 
dense  brittle  rock  is  shown  in  figure  la.  In  the  early  stage  of 
the  test,  the  volumetric  change  is  negative  (i.e.,  volume 
decrease);  which  is  mainly  attributable  to  the  elastic  behavior 
of  the  rock  but  also  reflects  the  closing  of  some  open  cracks. 
At  about  50  percent  of  the  peak  stress  the  curve  of  the  volu¬ 
metric  strain  versus  axial  stress  starts  to  deviate  from 
linearity;  the  deviation  becoming  progressively  greater  with 
increasing  stress.  Eventually,  at  fracture,  the  volumetric 
curve  shows  a  net  volvune  increase  with  respect  to  the  original 
unstressed  configuration.  In  figure  lb,  the  data  have  been  pre¬ 
sented  differently,  in  the  form  of  the  variation  of  the  volu¬ 
metric  strain  with  respect  to  the  axial  strain.  This  curve 
shows  that  before  the  peak  stress  is  reached  the  accumulated 
inelastic  volume  increase  (i.e.,  the  difference  between  total 
and  elastic  volume  increase)  remains  relatively  small. 

The  phenomenon  of  dilatancy  in  rocks,  which  finds  its 
cause  in  several  mechanisms  of  fracturation,  is  controlled  by 
the  mean  pressure  the  initial  porosity  of  the  rock.  In  brittle 
rock  the  dilatant  behavior  observed  during  compression  test  is 
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Figure  1.  Results  of  a  hypothetical  triaxial  test  on 
rock  sample. 


associated  with  microcracking  and  propagation  of  cracks  parallel 
to  the  direction  of  maximum  compressive  stress  (Brace  et  al., 
1966;  Cook,  1970).  In  porous  sedimentary  rocks  dilatancy  is  not 
only  caused  by  an  increase  of  the  crack  space  (porosity)  but 
also  by  sliding  along  intergranular  surfaces  (cataclastic  flow). 
At  high  confining  pressure,  tendency  of  dilatancy  is  suppressed 
or  even  reversed  (Swanson  and  Brown,  1972),  depending  on  the 
initial  porosity  of  the  rock,  as  the  high  pressure  forces  intra¬ 
crystalline  flow  to  occur.  For  very  porous  rock,  negative  dila¬ 
tancy  can  even  occur  at  failure,  as  two  processes  compete: 
(1)  dilatancy  caused  by  shearing,  (2)  volume  decrease  caused  by 
collapse  of  the  pore  structure. 

Dilatation  is  thus  an  all  pervasive  property  of  hard 
rock,  and  the  stress  values  at  which  it  starts  reflect  more 
permanent  changes  occuring  in  the  rock  structure.  It  has 
finally  to  be  noted  that,  although  dilatancy  is  closely  associ¬ 
ated  with  the  process  of  macrofracturation,  there  is  no  close 
correlation  -  as  it  is  sometimes  speculated  -  between  suppres¬ 
sion  of  dilatancy  and  the  transition  between  brittle  to  ductile 
failure  (Edmond  and  Paterson,  1972). 

2-2  CRITICAL  REVIEW  OF  LABORATORY  EXPERIMENTS. 

Unfortunately,  only  a  few  of  the  investigations  on  the 
dilatancy  of  rocks  can  be  used  in  the  understanding  of  the  rock 
response  around  deep  underground  excavations.  Indeed,  many 
experimental  studies  of  rock  dilatation  find  their  motivation  in 
the  analysis  of  earthquake  precursors,  and  are  thus  concerned 
with  measuring  rock  dilatancy  at  very  high  confining  pressures 
(e.g..  Brace  et  al .  ,  1966;  Edmond  and  Paterson,  1972;  Shock  et 
al . ,  1973).  In  contrast,  modeling  the  response  of  rock  tunnels 
requires  data  on  rock  dilatancy  at  relatively  low  confining 
pressure  (of  the  orders  of  tens  of  bars,  as  opposed  to  the 
thousands  of  bars  which  have  been  reached  in  high  confining 
pressure  experiments).  Also  rocks  around  excavations  experience 
a  stress  path  which  involves  increase  of  the  deviatoric  stress 
accompanying  unloading  of  the  mean  pressure,  while  most  labora¬ 
tory  experiments  are  characterized  by  a  concomitant  increase  of 


both  the  confining  pressure  and  the  deviatoric  stress.  Finally, 
few  data  are  available  beyond  peak  strength  because  many  experi¬ 
ments  were  carried  out  on  a  "soft"  testing  machine.  (Actually, 
not  only  a  knowledge  of  the  rock  volume  change  in  the  post¬ 
failure  stage  is  required,  but  also  in  test  conditions  where 
failure  is  pervasive  throughout  the  rock  sample  —  the  kinematic 
constraints  during  triaxial  test  experiments  allows  localized 
failure  modes  to  develop.) 

2-3  THEORETICAL  MODELS  FOR  ROCK  DILATANCY. 

The  modern  approach  for  modeling  the  response  of 
geomaterials  is  based  on  the  theories  of  incremental  elasto- 
plasticity  which  involve  the  existence  of  a  yield  function,  f, 
and  a  plastic  potential,  g.  The  yield  function,  f,  marks  the 
boundary  of  the  elastic  state  in  the  stress  space;  f  is  not  only 
a  function  of  the  stress  but  also  of  some  measure  of  the  plastic 
defoirmation,  either  the  plastic  work  or  the  accumulated  plastic 
shear  strain.  The  response  to  an  increment  of  stress  dt  is 
elastic  if  the  stress  point  x  is  inside  the  current  yield 
surface,  or  if  the  current  stress  point  is  on  the  yield  surface, 
and  the  stress  increment  is  pointing  inside  the  yield  surface 
(elastic  unloading).  During  continued  plastic  flow,  i.e., 
during  a  loading  history  where  the  stress  point  remains  on  the 
yield  surface,  the  strain  increment  de  associated  with  the 
stress  increment  dx  is  compounded  of  an  elastic  part  de  (related 
to  da  by  Hooke’s  law)  and  a  plastic  part  de^  which  is  propor¬ 
tional  to  the  gradient  to  the  potential  surface  (flow  rule). 
The  condition  of  continued  plastic  flow  is  expressed  mathemat¬ 
ically  by 


df  =  0  (1) 

while  the  flow  rule  is  given  by 

deP  =  \  (2) 

Like  the  yield  function,  f,  the  plastic  potential,  g,  is  a 
function  of  the  stress  and  of  some  measure  of  past  plastic  flow. 
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In  this  investigation,  we  are  dealing  with  material 
models  characterized  by  Mohr-Coulomb  yield  and  potential  func¬ 
tions  intended  to  simulate  the  response  of  pressure-sensitive 
dilatant  rock  materials.  The  Mohr-Coulomb  functions  are  of  the 
intrinsic  curve  type;  i.e.,  the  plastic  deformation  is  indepen¬ 
dent  of  the  intermediate  principal  stress  T2*  For  a  Mohr- 
Coulomb  material,  the  constitutive  equations  for  continuous 
plastic  flow  reduce  to 

ds  -  pdP  =  hdv  (3) 

dA  =  pdY  (4) 

where  dP  =  (dx^  +  dT2)/2,  ds  =  (dx^^  -  dx2)/2,  and  dA  and  dy 
represent,  respectively,  the  variation  of  plastic  volume  change 
(de^  +  de^)  and  plastic  distortion  (de^  -  de^)/  and  where  h 
represents  the  hardening  modulus.  Figure  2  gives  the  geo¬ 
metrical  interpretation  of  the  dilatancy  parameter  p  =  sin  (|)* 
and  the  friction  coefficient  m  =  sin<)).  Note  that  if  the  0  =  $*, 
the  flow  rule  is  associated  ($  represents  the  highest  theoreti¬ 
cal  value  for  the  dilatancy  angle 

In  the  previous  investigations  of  elastoplastic  models 
of  tunnels  (Labreche  and  Auld,  1980;  Reed  et  al . ,  1983;  AA, 
1983)  it  is  assumed  that  (1)  the  coefficient  of  friction  is  a 
constant,  (2)  the  material  is  nonhardening;  i.e.,  h  =  0,  and 
(3)  the  dilatancy  factor  is  a  constant.  These  assumptions  imply 
that  both  yield  and  potential  functions  are  fixed  in  the  stress 
space,  and  that  only  three  parameters  are  needed  to  describe  the 
plastic  deformation  (q,  the  unconfined  compressive  strength,  the 
friction  angle  <() ,  and  the  dilatancy  angle  0*).  However,  as 
discussed  in  the  introduction,  the  weakest  of  these  assumptions 
is  the  hypothesis  of  a  constant  dilatancy  angle,  since  it 
results  in  unbounded  inelastic  volume  changes.  In  the  present 
investigation,  the  first  two  assumptions  are  maintained  while 
the  third  assumption  is  relaxed  to  a  variable  condition. 

The  experimental  results  reported  in  the  previous 
section  suggest  that  the  dilatancy  parameter  is  both  a  function 


Figure  2. 
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of  the  accumulated  plastic  shear  strain  and  the  mean  pressure. 
However,  since  we  are  interested  in  the  phenomenon  of  dilatancy 
at  relatively  low  confining  pressure,  and  in  the  interest  of 
keeping  the  model  simple,  the  dependence  upon  the  mean  pressure 
will  here  be  ignored. 

The  proposed  flow  rule  is  based  on  an  exponential 
decay  of  the  dilatancy  factor  K*  (which  is  defined  as  minus  the 
ratio  of  the  maximum  to  the  minimum  plastic  strain  rates,  i.e., 

’'J '  ^  ^  ‘"Sj  •  “p  (■ 

This  law  is  based  on  some  limited  experimental  evi¬ 
dence  which  suggests  that  (1)  the  rate  of  dilatation  at  peak 
stress  in  dense  brittle  material  appears  to  be  consistent  with 
an  associated  flow  rule  (Ladanyi  and  Don,  1970;  Gerionnopoulos 
and  Brown,  1978),  and  that  (2)  the  rate  of  dilatation  progres¬ 
sively  drops  to  zero,  beyond  the  peak  stress.  The  parameters  y* 
in  the  exponential  law  (equation  5)  can  most  usefully  be  related 
to  the  maximum  inelastic  volume  increase  A*  by  integrating  the 
relation 


^  =  2e _ _ 

dy  KJ  +  1 

to  yield 

K  +  1 
=  V*  - 


(6) 


(7) 


Section  3  which  follows  outlines  the  development  of  an  elasto- 
plastic  model  of  a  deep  tunnel  which  is  based  on  the  variable 
flow  rule  (equation  5). 
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SECTION  3 

ELASTOPLASTIC  MODELS  OF  TUNNELS  UNDER 
HYDROSTATIC  AND  NONHYDROSTATIC  LOADING 

3-1  INTRODUCTION. 

In  this  chapter,  the  development  of  two  elastoplastic 
models  of  a  deep  cylindrical  tunnel,  which  implement  the 
material  model  described  in  Section  2  is  outlined.  Two  models 
are  considered,  one  for  hydrostatic  far-field  loading,  the  other 
one  for  nonhydrostatic  loading.  The  nonhydrostatic  model  is 
restricted  to  loading  conditions  for  which  the  problem  remains 
statically  determinate;  as  is  always  the  case  for  the  hydro¬ 
static  loading.  This  restriction  ensures  that  many  features  of 
the  solution  that  were  derived  for  a  constant  dilatancy  angle 
still  apply  for  the  improved  dilatation  model  (e.g.,  the  extent 
and  shape  of  the  failed  region  around  the  tunnel ) . 

The  mathematical  foundation  of  the  elastoplastic 
models  is  developed  for  a  hydrostatic  case  in  Appendix  A  and  in 
Appendixes  B  and  C  for  the  nonhydrostatic  far-field  loading. 
Appendix  B  details  a  formulation  for  calculating  the  tunnel 
closure  for  the  case  of  a  constant  dilatancy  angle  that  is 
developed  for  solving  the  general  case  with  variable  dilatancy 
in  Appendix  C.  The  formulation  discussed  in  Appendix  B  is  an 
alternative  to  that  derived  previously  (Detournay,  1983).  It 
was  developed  because  the  original  formulation  could  not  be 
implemented  easily  with  a  material  characterized  by  a  varicible 
dilatancy. 

The  equations  for  calculating  the  closure  of  the 
tunnel  are  derived  for  a  stress  history  intended  to  simulate  the 
excavation  unloading  of  a  prestressed  rock  mass.  As  discussed 
in  the  AA  report  (1983),  an  elastic  correction  has  to  be 
applied,  to  account  for  a  stress  history  with  a  far-field  stress 
surcharge . 
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3-2 


HYDROSTATIC  MODEL. 


3-2.1  Symmetry  Conditions. 

In  the  hydrostatic  model,  the  boundary  conditions 
consist  of  an  internal  pressure  p  inside  the  tunnel  and  a  mean 
pressure  at  infinity  (see  figure  4).  Because  of  the  symmetry 
of  the  boundary  conditions  and  the  geometry  for  this  problem, 
the  stress,  strain,  and  displacement  field  in  the  medium  depend 
only  on  the  radial  coordinate  r  of  the  cylindrical  coordinate 
system  which  has  its  origin  at  the  center  of  the  cavity. 

Provided  that 


p  <  — -  /  P°  -  — ^ - \  -  — ^ -  ( 8 ) 

P  K  +  1  r  K  -  1  )  K  -  1 
P  \  P  /  P 

(it  is  assumed  that  P^  >  2q),  the  tunnel  is  surrounded  by  a 
plastic  zone  of  external  radius  aR^-  Since  the  problem  is 
statically  determinate,  the  normalized  radius  of  the  elasto- 
plastic  interface  is  only  a  function  of  p,  P^,  and  the  yield 
parameters  q  and  Kp: 

l/(Kp-l) 

(9) 


^o  = 


P°  + 


K  +  1  p  + 


(Kp  - 


(Kp  -  1) 


Since  only  the  loading  processes  of  interest  (either  internal 
unloading  or  external  loading)  cause  a  mono tonic  increase  of  the 
radius  of  the  plastic  zone,  it  is  advantageous  to  use  R^  instead 
of  either  p  or  P°  as  a  kinematic  parameter.  In  other  words,  to 
determine  the  mechanical  fields  (stress,  strain,  and  displace¬ 
ment)  as  a  dual  function  of  the  coordinate  r  and  the  history 
parameter  R^. 

The  problem  being  statically  determinate,  the  stress 
field  in  the  medium  and  the  displacement  field  in  the  elastic 
region  r  >  aR  are  actually  independent  of  the  flow  rule. 
(Thei’*  expressions  can,  for  example,  be  found  in  Newmark,  1970, 
or  Hendron  and  Aiyer,  1971.)  To  calculate  the  closure  of  the 
tunnel  as  a  function  of  R  ,  determine  the  displacement  field 


W  /  y  / 


u(r,  R^)  in  the  plastic  zone  (a  <  r  <  aR^)-  However,  as  proven 
in  Appendix  A,  the  general  foirm  of  the  displacement  field 
u(r,  R^)  is  given  by 


u(r,  R^) 


i.e.,  the  displacement  is  only  a  function  of  the  ratio  r/aR^. 

3-2.2  Concept  of  the  Unit-Plane. 

Conceptually,  it  is  advantageous  to  introduce  the  unit 
plane  transformation  p  =  r/aR^.  In  the  unit  plane,  the  circle 
of  radius  p  =  1  separates  two  regions;  an  interior  plastic  one 
from  an  exterior  elastic  one.  As  the  plastic  annulus  grows  as  a 
consequence  of  changes  in  the  boundary  conditions,  the  image  of 
a  physical  point  in  the  unit  plane  moves  inwards  along  a  radial 
line,  crossing  the  unit  circle  when  the  radius  of  the  plastic 
boundary  reaches  that  physical  point. 

The  concept  of  the  unit  plane  is  a  powerful  one, 
because  it  substantially  simplifies  the  mathematics  of  the 
problem. 

3-2.3  Governing  Equations. 

The  normalized  displacement  u(p)  (p  <  1)  is  calculated 
by  solving  the  differential  equation 


•  vTF,  IT, . 
'w'.*  O  ‘A* 

lA.'AV.'S 


•  s  .s 


with 


p^u"  +  K*  pu'  -  K*  u  =  -A*p 


=  (K  -  1)  (K*  -  1)  +  (1  -  2v)  (K  +  1)  (K*  +  1) 


subject  to  the  boundary  conditions 

u(l)  =  -  1  ;  u'(l)  =  1  (12) 

The  dilatancy  factor  K*  is  a  function  of  the  maximum  plastic 
shear  distortion  -y  /  which  can  be  expressed  in  terms  of  p,  u,  and 
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The  differential  equation  can  thus  be  cast  in  the  general  form 

u"  =  F(p,  u,  u')  (14) 

which  is  well  suited  for  numerical  solution,  using  an  algorithm 
such  as  Runge-Kutta.  (Such  an  algorithm  is  usually  included  in 
the  math  library  of  a  scientific  calculator,  thus  making  the 
solution  of  equation  14  straightforward. ) 

The  function  u(p)  depends  only  on  three  dimensionless 
parameters:  Kp,  u ,  and  A*,  which  is  defined  as 

A*  =  ^  A*  (15) 

Experimental  evidence  suggests  that  the  maximum  inelastic  volume 
increase  A*  is  generally  less  than  5  percent,  thus  leading  to  a 
possible  range  of  values  for  A*  between  0  and  100. 

Once  the  function  u(p)  has  been  determined,  the  dis¬ 
placement  at  the  tunnel  wall,  as  a  function  of  the  history 
parameter  R^,  is  simply  given  by 

aR  S .  / 1  \ 

“  =  “  (y 

The  variation  of  the  normalized  displacement  (2G/aS°)  u  at  the 

boundary  has  been  plotted  in  figure  5  as  a  function  of  the 

normalized  radius  of  the  plastic  zone  for  K  =  3,  v  =  0.25,  and 

^  P 

various  values  of  the  parameter  A^.  To  illustrate  how  the  use 
of  a  constant  dilatancy  angle  can  be  misleading,  figure  6  pre¬ 
sents  the  apparent  constant  dilatancy  factor  K*  which  yields, 
for  any  given  value  of  R^,  the  same  displacement  at  the  boundary 
as  obtained  with  the  vari€d)le  flow  rule  (equation  5).  Figure  7 
suggests  that  the  solution  for  A^=  0.1  should  approximate 
closely  the  solution  for  a  plastically  incompressible  material 
(K*  =  1)  for  values  of  R^  greater  than  1.5;  while  the  displace¬ 
ment  for  A.=  100  should  be  close  to  the  one  predicted  by  an 
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associated  flow  rule  (K*  =  3)  for  R  less  than  2.5.  The  numeri- 

p  o 

cal  analysis  indicates  that  the  plastic  dilatation  on  the  boun¬ 
dary  rapidly  reaches  its  maximum  value  in  the  first  case 
(A  s  0.1  for  R  =1.1;  A*  =  0.1),  but  that  only  45  percent  of 
the  maximvun  inelastic  dilatation  has  been  reached  for  R^  =  2.5 
in  the  second  case  (A*  =  100). 

3-3  NONHYDROSTATIC  MODEL. 

3-3.1  Modes  of  Failure. 

In  this  model  the  stress  at  infinity  is  characterized 
by  a  mean  pressure  P°  and  a  stress  deviatoric  S°  (see  figure  7). 
In  this  case,  the  problem  is  characterized  by  two  axes  of  sym¬ 
metry,  which  are  parallel  to  the  principal  stress  directions  at 
infinity.  Due  to  the  existence  of  a  stress  deviatoric  S°  at 
infinity,  different  modes  of  failures  can  develop  around  the 
tunnel  depending  on  the  relative  values  of  P°,  S°,  p,  and  q. 

Consider  first  the  case  of  an  unsupported  tunnel  (p  =  0).  The 
different  types  of  behaviors  can  graphically  be  depicted  in  the 
normalized  stress  space  (P°/q,  S*^/q)  (see  figure  8); 

•  Type  1:  Elastic  behavior  only  (region  desig¬ 

nated  I ) 

•  Type  2 ;  Limited  failure  in  a  direction  perpen¬ 

dicular  to  the  major  in-situ  stress 
(region  designated  I la) 

•  Type  3 :  Tunnel  completely  surounded  by  an 

oval-shaped  yield  zone  (region  desig¬ 
nated  Ilb) 

•  Type  4:  A  "butterfly" -shaped  plastic  region 

around  the  tunnel  (region  designated 
III) 

Region  II  in  the  stress  diagram  corresponds  to  statically  deter¬ 
minate  cases;  i.e.,  conditions  for  which  the  extent  and  shape  of 
the  plastic  region  are  entirely  controlled  by  the  stress 
boundary  conditions  and  the  yield  parameters  q  and  K^.  The 
boundary  between  statically  determinate  and  indeterminate  condi¬ 
tions  corresponds  to  a  line  of  critical  obliquity  m*  (recall 
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that  the  obliquity  is  defined  as  the  ratio  of  S°/S°),  The 
critical  obliquity  m*  is  only  a  function  of  the  friction  angle  0 
(see  table  1).  Far-field  stress  points  on  a  line  of  equal 
obliquity  correspond  to  elastoplastic  interfaces  that  are  geo¬ 
metrically  identical,  but  of  different  sizes. 


TABLE  1.  Critical  obliquity  m*. 


4> 

0“ 

(-• 

0 

0 

0 

0 

OJ 

30“ 

0 

0 

0.414 

0.437 

0.466 

0.500 

0.542 

The  existence  of  an  internal  support  pressure  in  the 
tunnel  causes  the  boundary  between  regions  I  and  I la  and  the 
boundary  between  regions  I  la  and  lib  to  move  to  the  left.  As 
discussed  in  AA  (1983),  the  effect  of  an  internal  support  pres¬ 
sure  can  be  taken  into  account  by  a  simple  geometrical  construc¬ 
tion,  which  involves  moving  the  far-field  stress  point  along  a 
line  of  equal  obliquity.  Accordingly,  the  presence  of  an 
internal  pressure  changes  neither  the  boundary  between  stati¬ 
cally  determinate  and  indeterminate  conditions,  nor  does  it 
change  the  shape  of  the  elastoplastic  interface. 

3-3.2  Consequences  of  Statical  Determinacy. 

The  semianalytical  solution  developed  in  this  report 
is  restricted  to  far-field  conditions  for  which  the  problem  is 
statically  determinate;  i.e.,  for  obliquity  less  than  the  criti¬ 
cal  obliquity  ra* .  The  restriction  to  statically  determinate 
conditions  ensures  that,  similarly  to  the  hydrostatic  loading, 
the  location  of  the  elastoplastic  boundary  is  solely  controlled 
by  the  stress  boundary  conditions  P°,  S°,  and  p,  and  the  yield 
parameters  q  and  Kp.  For  conditions  for  which  the  tunnel  is 
completely  engulfed  by  a  plastic  zone,  the  interface  is  char¬ 
acterized  by  a  major  to  minor  axis  ratio  equal  to 


major  axis  _  / 1  +  m  ^ 

minor  axis  \1  -  m/ 

with  the  obliquity  m  defined  as 


(17) 


(18) 


The  average  radius  of  the  plastic  zone  corresponds,  however,  to 
the  one  computed  for  hydrostatic  loading.  Also,  the  stress 
field  in  the  medium  and  the  displacement  field  in  the  elastic 
region  are  the  same  as  the  one  computed  for  a  constant  dilatancy 
angle.  To  calculate  closure  of  the  tunnel,  we  need  to  calculate 
the  displacement  in  the  plastic  region  as  a  function  of  the 
boundary  conditions,  or  equivalently  as  a  function  of  the 
average  radius  of  the  plastic  zone  (as  for  the  case  of  hydro¬ 
static  loading).  In  this  case,  however,  the  displacement  has 
not  only  a  radial  component  but  also  a  tangential  one  (except  on 
the  axes  of  symmetry),  and  is  a  function  of  both  cylindrical 
coordinate’s  (r,0)  and  of  the  kinematic  parameter  R^. 

3-3.3  Governing  equations. 

Derivation  of  the  equations  governing  the  displacement 
in  the  plastic  region  of  a  material  with  a  variable  dilatancy 
follows  closely  the  approach  adopted  for  a  constant  dilatancy 
( see  Appendixes  B  and  C ) .  The  derivation  is  based  on  the  pro¬ 
perty  exhibited  by  the  stress  solution  that  there  is  no  rotation 
of  the  principal  stress  directions  in  the  plastic  zone  during 
propagation  of  the  failure  zones  around  the  tunnel .  This  pro¬ 
perty  of  the  solution  ensures  that  the  principal  directions  of 
the  incremental  plastic  strain  tensor  de^  remains  locked  in  the 
radial  and  tangential  directions  everywhere  in  the  plastic  zone 
and  at  any  time  during  the  monotonic  loading.  This  feature  of 
the  problem  allows  integration  of  the  flow  rule 


(19) 


2  +  (Kp  -  l)e  "  * 


with 


Y*  = 


Figure  9  illustrates  the  geometrical  interpretation  of  the 
tangent  and  secant  dilatancy  factor  (note  that  for  a  constant 
dilatancy,  the  distinction  between  tangent  and  secant  dilatancy 
factor  disappears ) . 

The  governing  equations  of  the  displacement  in  the 
plastic  region  are  deduced  from  equation  (20)  as  follows:  the 
plastic  strain  is  expressed  as  the  difference  between  the  total 
and  the  elastic  strain;  the  elastic  strain  field  in  the  plastic 
zone  is  explicitly  determined,  using  Hooke's  law  and  the  stress 
solution;  and,  the  total  strain  can  be  related  to  the  partial 
spatial  derivatives  of  the  displacement.  These  relationships 
result  in  a  set  of  partial  differential  equations  for  the  dis¬ 
placement  that  are  solved  using  the  elastic  displacement  on  the 
elastoplastic  interface  as  boundary  conditions.  It  turns  out, 
that  as  for  the  hydrostatic  case,  the  general  form  of  the  dis¬ 
placement  is  given  by 


U  (r,  6;  R^)  = 


SI  . 


The  displacement  field  is  thus  a  unique  function  of  the  coordi¬ 
nates  of  the  unit-plane. 

All  calculations  done,  the  partial  differential  equa¬ 
tions  governing  the  normalized  displacement  in  the  plastic  zone 
are  given  by 

2*  -  sin  i.  =  H,  (p.t) 

sin  24.  =  (P.P) 

'  '  /  A  V 


where 


with 


2X*  K  -1 

(p,<|>)  =  -  -  p  cos  2(1) 

(Xp  -  1)  (XJ  -  1) 

X  +  1  X»  -  1 

+  2(l-2\>)  7^ - r  cos  20  +  2m  - 

S  ^  KJ  +  1 

^2  ~  2m  sin  20 


=  (Kp  -  1)(KJ  -  1)  +  (1  -  2v)(Kp  +  1)(K*  +  1)  (25) 


This  nonlinear  system  of  equations  is  of  the  hyperbolic  type, 
and  can  thus  be  solved  by  the  method  of  characteristics.  This 
system  is  however  very  stiff  and  as  a  result,  poor  accuracy  is 
achieved  if  standard  algorithms  of  solution  are  used.  Because 
of  the  ill-conditioned  nature  of  equation  (23),  a  special  scheme 
had  to  be  used.  This  involved  tranforming  the  equations  into  a 
system  of  two  ordinary  differential  equations  along  the  char¬ 
acteristics  and  solving  them  using  a  central  node  finite  dif¬ 
ference  technicjue. 
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3-3.4  Validation  and  Application. 

As  a  means  of  validating  the  numerical  algorithm 
(program  KINVAR)  developed  to  calculate  the  closure  of  a  tunnel 
in  the  presence  of  a  nonhydrostatic  far-field  stress,  two  hydro¬ 
static  test  cases  were  run  and  compared  with  results  obtained  by 
the  code  GROUND,  which  solves  the  differential  equation  with  the 
Runge-Kutta  algorithm.  The  comparison  is  shown  in  figures  10 
and  11,  where  the  normalized  radial  displacement  at  the  tunnel 
wall  is  plotted  as  a  function  of  the  radius  R^  of  the  elasto- 
plastic  interface.  These  figures  indicate  that  with  relatively 
few  points,  (which  is  indicative  of  a  relatively  coarse 
characteristic  mesh),  KINVAR  is  able  to  faithfully  predict  the 
closure  of  the  tunnel . 

Two  nonhydrostatic  cases  have  also  been  solved,  and 
the  results  shown  in  figures  12  and  13  (radial  displacement  on 
the  two  axes  of  symmetry  as  a  function  of  These  plots 
confirm  that  the  tunnel  ovals  during  closure,  with  the  direction 
of  maximum  closure  becoming  perpendicular  to  the  maximum  com¬ 
pressive  far-field  stress,  once  the  plastic  region  around  the 
tunnel  becomes  sufficiently  large. 
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SECTION  4 

LABORATORY  EXPERIMENTS 

4-1  MODEL  STUDIES. 

Over  the  past  15  years  or  more  the  Department  of 
Defense,  through  several  agencies  and  contractors,  has  conducted 
a  large  number  of  scaled  model  studies  of  underground  excava¬ 
tions,  mostly  in  rock  simulant.  The  studies  have  investigated  a 
number  of  significant  parameters  governing  the  behavior  of 
hardened  excavations,  including:  alternative  ground  support 
systems,  monotonic  loading,  and  static  versus  dynamic  loading. 
The  tests  have  undoubtedly  contributed  significantly  to  the 
understanding  of  the  behavior  of  underground  excavations  and  the 
results  of  the  scaled  model  tests  have  been  demonstrated  to  be 
qualitatively  similar  to  those  observed  during  testing  of  full- 
scale  structures  in  the  vicinity  of  underground  weapons  tests . 

The  most  extensive  series  of  tests  have  been  performed 
by  the  Stanford  Research  Institute  (SRI),  using  cylindrical 
specimens  4  in.  and  12  in.  in  diameter  into  which  excavations 
have  been  drilled  or  cast.  Emphasis  during  the  earlier  tests 
was  on  simulating  the  behavior  of  excavations  in  tuff,  and  a 
niunber  of  rock  simulants  with  relatively  low  friction  angles 
were  used  to  fabricate  lined  and  unlined  tunnels  in  intact  and 
jointed  rock.  Emphasis  was  also  placed  on  understanding  the 
impact  of  static  versus  dynamic  loading.  More  recently  the 
interest  in  siting  a  deep  based  missile  system  in  a  sedimentary 
stratigraphy  has  resulted  in  testing  of  a  number  of  different 
types  of  structures  in  rock  simulants  exhibiting  higher  friction 
angles . 

As  part  of  the  present  investigation  a  detailed  review 
of  the  scaled  model  tests  was  conducted,  with  a  view  to  identi¬ 
fying  data  that  may  be  used  to  validate  the  variable  dilatancy 
model  described  in  the  previous  sections  of  this  report. 
Despite  the  large  number  of  tests  that  have  been  performed,  very 
few  lend  themselves  easily  to  the  desired  purpose.  There  are 
two  important  reasons  for  this  difficulty.  The  first  is  that 
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the  majority  of  the  tests  have  been  performed  under  conditions 
intended  to  reproduce  the  uniaxial  strain  loading  that  might  be 
experienced  by  an  excavation  as  a  consequence  of  a  distant 
nuclear  burst,  for  example.  As  will  be  discussed  later,  the 
selected  loading  path  rapidly  leads  to  conditions  that  are 
statically  indeterminate  and  are,  therefore,  not  amenable  to 
analysis  in  closed  form.  The  second  problem  is  that  the  region 
of  failed  rock  around  the  tunnels  is  typically  relatively  large 
when  compared  to  the  dimensions  of  the  cylinder  within  which  it 
is  located.  Under  such  circumstances  the  influence  of  the 
boundary  can  significantly  modify  the  response.  Again,  this 
effect  cannot  be  accounted  for  in  the  analytical  model .  The 
third  problem  is  that  most  of  the  support  systems  tested  exert  a 
nonuniform  pressure  on  the  rock  simulant  unless  the  model  is 
subjected  to  isotropic  loading.  As  currently  developed,  the 
analytical  model  requires  that  the  support  be  approximated  as  a 
uniform  internal  pressure.  The  approximation  may  be  acceptable 
for  rock  bolted  or  backpacked  structures,  but  is  inappropriate 
for  integral  steel  or  concrete  liners  subjected  to  nonisotropic 
loading. 

In  the  following  sections  the  results  of  a  number  of 
laboratory  tests  are  discussed.  Considering  the  problems 
associated  with  the  uniaxial  loading  conditions,  emphasis  is 
placed  on  cases  in  which  the  loading  was  isotropic.  However, 
there  is  also  a  discussion  of  the  results  of  selected  tests 
under  uniaxial  strain  conditions. 

4-2  ISOTROPIC  LOADING  OF  LOW  FRICTION  SIMULANTS. 

A  large  number  of  tests  have  been  performed  using  a 
tuff  simulant  designated  RMG-2C2 .  Typical  results  from  early 
tests  on  lined  and  unlined  tunnels  in  intact  rock  subjected  to 
static  and  dynamic  loads  are  reproduced  in  figure  14.  These 
tests  revealed  a  significant  difference  between  behavior  under 
static  and  dynamic  loads,  and  dry  and  saturated  conditions.  The 
differences  were  attributed,  in  the  most  part,  to  pore  water 
effects;  with  the  pore  water  weakening  the  specimens  in  the 
static  tests  and  strengthening  them  in  the  dynamic  tests . 
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Following  this  finding,  testing  has  emphasized  static  loading  of 
saturated  samples.  Under  such  circumstances  the  properties  of 
the  RMG-2C2  simulant  are  reported  to  be  (SRI,  1979): 

Young's  Modulus  E  =  1.6*10  psi 

Poisson's  Ratio  u  =  0.18 

Unconfined  Compressive  Strength  q  =  3200  psi 
Friction  Angle  $  =  2.5° 

For  reasons  discussed  above  we  shall  consider  the  case  of  lined 

excavations  subjected  to  isotropic  loading.  Typical  results  of 
such  tests  conducted  with  different  thicknesses  of  A1  6061 
aluminum  liner  are  reproduced  in  figure  15.  Superimposed  on  the 
laboratory  data  are  closure  versus  applied  pressure  curves 
computed  using  a  model  that  assumes  associated  behavior  (4»  = 
for  the  simulant  and  the  strain  hardening  property  illustrated 
in  figure  16.  It  may  be  observed  that  the  analytical  model 
reasonably  reproduces  the  laboratory  data.  However,  as  dis¬ 
cussed  below,  this  finding  cannot  be  regarded  as  evidence  of  the 
validity  of  the  constitutive  model  assumed  for  the  rock 
simulant. 

The  most  important  factor  affecting  the  value  of  this 
data  is  the  very  low  friction  angle.  In  figure  17  the  tunnel 
closure  histories  predicted  using  full  dilatation  ((|)*  =  0)  and 
zero  dilatation  (0*  =  0)  are  illustrated  for  lined  and  unlined 
tunnels  in  the  RMG-2C2  simulant.  (Models  based  on  these  two 
extreme  assumptions  have  been  referred  to  as  the  Hendron  and 
Newmark  models  respectively.)  If  figure  17  is  compared  to 
figure  15  it  is  clear  that  the  differentiation  between  the  full 
dilatation  and  zero  dilatation  is  not  an  important  phenomenon 
for  such  a  low  friction  material . 

The  second  consideration  for  these  tests  on  low  fric¬ 
tion  materials  is  the  boundary  conditions.  For  the  isotropic 
loading  conditions  considered  here  this  effect  can  be  investi¬ 
gated  analytically  using  the  solution  for  a  thick-walled 
cylinder  subjected  to  uniform  internal  and  external  pressures  (p 
and  Pq).  The  extent  of  the  plastic  region  in  the  cylinder  is 
given  by  the  equation  (Kennedy,  1975): 
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in  which  a  and  b  are  respectively  the  internal  and  external 
radii  of  the  cylinder  and  R  is  the  radius  of  the  plastic  zone 
normalized  by  a.  As  noted  earlier,  the  radius  of  the  plastic 
zone  around  a  tunnel  in  an  Infinite  region  is  given  by: 


R  = 
o 


P°  + 


2 

K_  +  1 


(K^  -  1) 


(Kp  -  1) 


l/(Kp-l) 


which  is  a  special  case  of  the  thick  walled  cylinder  equation, 
that  can  be  deduced  by  allowing  (b/a)  to  approach  infinity. 
Equation  (26)  can  be  solved  to  define  the  minimum  thickness 
below  which  the  thick-walled  cylinder  will  be  completely  plas¬ 
tic.  This  thickness  is  a  function  of  the  internal  and  external 
pressures  and  the  material  properties: 


(i).  = 

'  'min 


p°  + 


(K  -  1) 


(Kp  -  1) 


l/(Kp-l) 


Referring  to  equation  (27),  it  can  be  seen  that  this 
minimum  thickness  can  be  much  larger  than  the  extent  (R^)  of  the 
plastic  zone  around  a  hole  in  an  infinite  region.  Specifically: 


(“lin  l) 


l/<Kp-l) 


friction 


Clearly  the  difference  is  more  important  for  very  low 
angle  materials,  such  as  the  RMG-2C2  rock  simulant. 
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This  can  be  illustrated  best  by  solving  equation  (26)  to  obtain 
the  relationship  between  the  radius  of  the  plastic  zone  and  the 
dimensions  of  the  thick-walled  cylinder  for  a  particular  choice 
of  material  properties.  Such  a  relationship  is  illustrated  in 
figure  18.  That  figure  is  for  the  case  of  a  thick-walled  cylin¬ 
der  of  RMG-2C2  without  internal  pressure  and  with  external 
pressure  expressed  as  a  multiple  of  the  uniaxial  strength.  At 
relatively  high  loads  the  extent  of  the  plastic  zone  in  the 
finite  region  is  considereibly  enlarged,  which  indicates  that  the 
displacements  observed  in  these  models  will  be  much  greater  than 
would  be  observed  in  the  field.  Hence  the  results  of  tests  of 
tunnels  in  RMG-2C2  simulant  are  unsatisfactory  for  validating 
the  variable  dilatancy  model  on  two  counts.  First,  the  friction 
angle  is  too  small  to  provide  a  differentiation  between  alterna¬ 
tive  dilatation  models.  Second,  the  results  of  the  model 
studies  are  strongly  influenced  by  the  boundary  conditions. 

4-3  ISOTROPIC  LOADING  OF  HIGH  FRICTION  SIMULANTS. 

A  series  of  tests  on  lined  and  unlined  tunnels  were 
performed  by  SRI  on  a  relatively  high  friction  simulant  desig¬ 
nated  6B.  Material  properties  of  this  simulant  are  reported 
(Lindberg,  1983): 

Young's  Modulus  E  =  2.0*10  psi 

Poisson's  Ratio  v  =  0.25 

Unconfined  Compressive  Strength  q  =  4300.  psi 

Friction  Angle  (J)  =  33° 

Results  for  four  cases  of  isotropic  loading  are  reproduced  in 
figure  19  (Lindberg,  1983).  Superimposed  on  the  experimental 
data  are  theoretical  results  obtained  using  the  full  (associ¬ 
ated)  dilatation  model.  The  theoretical  results  are  for  three 
different  steel  liner  thicknesses  in  addition  to  the  unlined 
case . 

It  may  be  observed  from  figure  19  that  there  are 
significant  differences  between  the  experimental  data  and  the 
theoretical  model.  The  most  obvious  is  that  the  effect  of  the 
liner  is  overestimated  in  all  cases. 


OUTER  RADIUS  OF  CYLINOER/HOLE  RADIUS,  b/a 


Figure  18.  Relationship  between  radius  of  plastic  region  around 
hole  and  radius  of  thick-walled  cylinder.  Curves 
are  plotted  for  selected  normalized  external  pressures 
for  material  with  a  friction  angle  of  2.5°. 
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This  is  relatively  easy  to  explain  for  the  very  thin  liner  since 
it  was  reported  that  the  thin  liners  were  severely  buckled 
during  loading.  For  the  thicker  liners  it  must  be  assumed  that 
their  effect  is  delayed  by  a  certain  amount  of  "consolidation” 
and  elastic  response  before  the  full  support  pressure  is 
mobilized,  when  the  liner  is  fully  yielding. 

A  second  important  departure  of  the  experimental 
behavior  from  the  theoretical  model  is  that  the  theoretical 
model  appears  to  overestimate  the  closure  at  higher  loads . 
Unfortunately  the  uncertainty  as  to  the  efficiency  of  the  liners 
makes  it  difficult  to  quantify  this  effect.  Despite  this  uncer¬ 
tainty  it  is  instructive  to  compare  these  test  results  with 
predictions  made  using  the  variable  dilatancy  model.  First, 
however,  it  is  appropriate  to  evaluate  the  importance  of  the 
boundary  conditions,  and  to  ascertain  whether  these  are  an 
important  consideration  in  this  case. 

In  figure  20  the  relationship  between  the  radius  of 
the  plastic  region  and  the  radius  of  a  thick-walled  cylinder  of 
6B  rock  simulation  subjected  to  external  pressure  is 
illustrated.  The  format  of  that  figure  is  the  same  as  fig¬ 
ure  18,  except  that  the  vertical  scale  has  been  enlarged  because 
the  effect  of  using  a  higher  friction  simulant  is  to  restrict 
the  growth  of  the  plastic  region.  From  the  figure  it  is  clear 
that  once  the  boundaries  lie  beyond  approximately  six  tunnel 
radii  they  cease  to  have  a  significant  influence.  Since  the 
tests  performed  by  SRI  satisfy  this  constraint  we  conclude  that, 
for  high  friction  simulants,  the  tunnel  deformation  should  be 
relatively  unaffected  by  the  fact  that  the  test  specimen  is 
finite . 

Figure  21  illustrates  the  relationship  between  tunnel 
closure  and  support  pressure  for  alternative  assumptions  regard¬ 
ing  the  maximum  inelastic  strain.  Superimposed  on  the  plots  are 
the  results  extracted  from  figure  19,  using  the  theoretical 
level  of  support  offered  by  the  three  liners.  It  is  evident 
that  either  the  effect  of  the  liners  is  being  overestimated  or 
that  there  is  an  initial  displacement  that  is  not  accounted  for 
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in  the  theoretical  model.  It  is  difficult  to  reconcile  either 
of  these  alternatives  with  the  fact  that  the  closures  appear  to 
be  underestimated  at  low  loads.  Despite  this  problem,  the 
results  of  this  comparison  do  tend  to  suggest  that  the  variable 
dilatancy  model  will  provide  an  improved  simulation  of  observed 
behavior.  There  is  evidence  that  the  rate  of  increase  of 
closure  does  decrease  at  the  higher  loads.  If  forced  to  esti¬ 
mate  the  maximum  inelastic  strain  using  this  limited  data  base 
it  might  be  set  at  approximately  10  percent. 

4-4  BIAXIAL  LOADING  OF  HIGH  FRICTION  SIMULANTS. 

Additional  data  for  the  response  of  lined  tunnels  in 
SRI  rock  simulant  6B  is  provided  by  Lindberg  (1983),  but  for  the 
case  of  simulated  uniaxial  loading.  In  this  case  the  solution 
for  the  closure  is  not  available  in  closed  form  because  the 
problem  rapidly  becomes  statically  indeterminate.  However,  in 
an  attempt  to  gain  further  data  on  the  properties  of  the  6B,  it 
was  considered  appropriate  to  perform  analysis  of  the  uniaxial 
loading  tests  using  a  finite  element  code  capable  of  simulating 
a  Mohr-Coulomb  material,  providing  associated  behavior  is 
assumed. 

Results  from  the  uniaxial  strain  model  studies  per¬ 
formed  by  SRI  are  reproduced  in  figure  22.  Before  attempting  to 
reproduce  those  results  a  preliminary  calculation  was  performed 
for  the  case  of  isotropic  loading.  The  result  of  that  calcula¬ 
tion  is  illustrated  in  figure  23,  where  it  is  used  to  verify 
that  the  numerical  prediction  was  in  good  agreement  with  the 
earlier  closed-form  solution.  Having  thus  confirmed  the  ade¬ 
quacy  of  the  numerical  model,  two  alternative  simulations  of  the 
uniaxial  strain  condition  were  investigated.  First,  uniaxial 
strain  was  imposed  by  restraining  lateral  displacement  of  the 
model.  To  be  reasoneibly  consistent  with  the  laboratory  configu¬ 
ration,  the  lateral  boundary  was  placed  6.5  tunnel  radii  from 
the  centerline  of  the  tunnel.  Since  this  boundary  is  compara¬ 
tively  close  to  the  tunnel,  it  is  probable  that  the  first  boun¬ 
dary  conditions  unrealistically  restrains  the  model.  Second,  a 
confining  stress  equal  to  that  generated  in  the  free  field  under 
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conditions  of  uniaxial  strain  was  imposed.  This  second  condi¬ 
tion  is  probably  too  compliant,  which  suggests  that  displacement 
of  the  tunnel  wall  would  be  overestimated. 

Results  from  the  two  numerical  simulations  of  uniaxial 
strain  loading  are  illustrated  in  figure  23.  Two  aspects  of  the 
behavior  are  very  interesting.  First,  the  results  are  extremely 
sensitive  to  the  assumed  boundary  conditions,  with  larger  dis¬ 
placements  resulting  from  the  stress  controlled  boundary. 
Second,  the  displacements  are  significantly  smaller  than 
observed  in  the  laboratory,  even  though  the  liners  were  not 
incorporated  in  the  numerical  simulation.  Since  the  stress 
controlled  boundary  conditions  should  have  resulted  in  an  over¬ 
estimate  of  the  tunnel  closure  it  appears  that  there  must  be 
some  deficiency  in  the  material  model.  Probably,  either  the 
reported  properties  are  incorrect,  or  the  Mohr-Coulomb  model 
does  not  adequately  describe  the  material  behavior.  These 
observations  are  reinforced  by  the  fact  that  the  numerical  model 
assumed  associated  behavior,  and  therefore  predicts  the  maximum 
possible  dilatation. 

4-5  EFFECT  OF  TUNNEL  REINFORCEMENT. 

SRI  provided  data  on  a  series  of  tests  to  evaluate  the 
effect  of  rockbolts  as  a  means  of  tunnel  hardening.  These  tests 
were  intended  to  simulate  18  ft  diameter  tunnels  either  unsup¬ 
ported  or  supported  with  #20  rockbolts  on  2  ft  centers.  This 
degree  of  reinforcement  amounts  to  an  internal  pressure  of 
approximately  680  psi  (4.7  MPa),  if  it  can  be  assumed  that  the 
bolts  exert  a  pressure  equal  to  the  yield  stress  of  the  steel . 
The  tests  selected  for  analysis  here  were  designated  by  SRI  as 
LSUX-35  and  LSUX-39.  The  reported  properties  of  the  rock  simu¬ 
lant  HF5  used  for  these  two  tests  are,  for  material  from  the  mix 
used  in  LSUX-39: 

Young's  Modulus  E  =  1.4*10  psi 

Poisson's  Ratio  u  =  0.25 

Unconfined  Compressive  Strength  q  =  4900  psi 

Friction  Angle  (j)  =  40° 
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Since  the  results  of  the  simulation  of  the  uniaxial 
loading  tests  performed  on  material  6B  indicated  a  considerable 
sensitivity  to  the  boundary  conditions,  the  load  path  used  by 
SRI  was  followed  as  closely  as  possible.  The  record  of  the 
relationship  between  the  vertical  and  lateral  pressures  during 
the  LSUX-39  test  is  reproduced  in  figure  24.  (This  relationship 
is  generated  by  slaving  the  lateral  pressure  to  maintain  zero 
circumferential  strain  at  sample  points  close  to  midheight  of 
the  cylinder  of  rock  simulant. )  To  simplify  the  numerical 
modeling,  this  load  path  was  idealized  as  three  linear  segments 
and  the  internal  support  pressure  was  applied  incrementally 
during  the  first  load  segment. 

Results  of  the  laboratory  tests  are  reproduced  in 
figure  25  and  those  of  the  numerical  simulation  in  figure  26 . 
In  both  cases  the  predicted  displacements  are  substantially  less 
than  observed  in  the  laboratory.  Since  any  uncertainty  in  the 
loading  conditions  was  removed  by  carefully  following  the 
laboratory  procedure,  it  is  concluded  that  there  are  deficien¬ 
cies  in  the  material  description.  Once  again,  this  may  be  in 
the  definition  of  the  properties  or  in  the  constitutive  model . 
The  most  likely  explanations  are  either  that  the  Mohr-Coulomb 
model  substantially  under-predicts  the  extent  of  the  plastic 
region  or  that  some  other  failure  mechanism,  such  as  near  sur¬ 
face  spalling,  is  occurring. 

4-6  CONCLUSIONS. 

The  results  of  evaluation  of  the  results  of  laboratory 
tests  using  models  based  on  the  closed-form  solution  and  a 
finite  element  procedure  indicate  that  the  Mohr-Coulotnb  model 
substantially  under-predicts  the  closure  even  when  associated 
behavior  is  assumed.  This  implies  that  the  material  description 
used  is  inadequate,  either  in  the  choice  of  material  properties 
or  in  the  constitutive  model .  Not  enough  is  known  about  the 
material  properties  (uniaxial  strength,  elastic  modulus, 
Poisson's  ratio,  and  friction  angle)  to  determine  whether  that 
is  the  source  of  error.  However,  it  is  reasonable  to  question 
whether  appropriate  account  has  been  taken  of  scale  effects  that 
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are  important  in  reality  but  are  ignored  in  the  mathematical 
descriptions . 

Given  the  uncertainty  in  the  material  properties  it  is 
difficult  to  draw  any  firm  conclusions  about  the  adequacy  of  the 
Mohr-Coulomb  model.  However,  the  fact  that  the  models  based  on 
associated  behavior  consistently  under  predict  the  observed 
deformation,  except  in  the  case  of  isotropic  loading,  suggests 
that  the  simple  constant  friction  plasticity  model  is  inade¬ 
quate.  Additional  laboratory  testing  would  be  required  to 
identify  the  nature  of  the  deficiency,  but  it  seems  most  likely 
that  the  model  underpredicts  the  extent  of  the  plastic  region 
around  the  tunnel.  Careful  inspection  of  cross  sections  of  a 
tunnel  structure  after  testing  could  be  used  as  a  means  to  test 
this  hypothesis.  Also,  consideration  could  be  given  to  moni¬ 
toring  microseismic  emissions  during  tests  to  detect  regions  of 
inelastic  behavior. 


SECTION  5 


CONCLUSIONS 

This  study  had  as  its  main  objective  the  development 
of  a  mathematical  model  of  a  deep  tunnel  based  on  an  improved 
dilatation  model  of  the  rock,  so  as  to  overcome  the  principal 
shortcoming  of  existing  models  based  on  the  assumption  of  a 
constant  dilatancy  angle. 

A  review  of  published  laboratory  experiments  on  rock 
dilatancy  revealed  that  very  few  experimental  data  could  be  used 
to  support  the  development  of  an  improved  dilatation  model  for 
rock  at  low  confining  pressure.  Reasons  for  this  deficiency  can 
be  found  in  the  use  of  "soft"  testing  machines,  which  are 
responsible  for  the  lack  of  data  beyond  the  peak  stress,  or  in 
the  fact  that  many  experiments  were  conducted  at  very  high 
confining  pressure,  so  as  to  simulate  the  behaviors  of  rock  at 
great  depth.  Some  limited  experimental  evidence  suggest,  how¬ 
ever,  that  at  peak  strength,  the  maximiun  theoretical  dilation 
rate  is  achieved  (associated  flow  rule)  and  that  it  progres¬ 
sively  diminishes  afterwards  with  the  plastic  shear  strain.  On 
that  basis,  a  very  simple  dilatation  model  was  implemented  which 
involved  the  introduction  of  a  single  parameter,  the  maximum 
inelastic  volume  increase,  in  contrast  to  the  constant  dilatancy 
angle  parameter  used  in  previous  dilatation  models.  This  con¬ 
stitutive  model  fits  between  the  two  extremes:  constant  dilata¬ 
tion  models  —  the  so-called  full  dilatation  model  (dilatancy 
angle  constant  and  equal  to  the  friction  angle)  —  and  the  zero 
dilatation  model  (zero  dilatancy  angle).  It  has  the  advantage 
of  relying  on  a  physically  meaningful  parameter. 

The  improved  dilatation  model  of  rock  was  then  used 
for  the  development  of  two  mathematical  models  of  a  deep  cylin¬ 
drical  tunnel,  one  for  hydrostatic,  the  other  one  for  nonhydro¬ 
static  loading.  For  the  hydrostatic  loading,  it  was  shown  that 
closure  of  the  tunnel  requires  the  solution  of  a  nonlinear 
ordinary  differential  equation,  and  for  the  nonhydrostatic 
loading  a  system  of  nonlinear  partial  differential  equations  of 


the  hyperbolic  type.  In  both  cases,  the  numerical  procedures 
are  discussed  in  detail:  Runge-Kutta  for  the  hydrostatic  load¬ 
ing,  and  the  method  of  characteristics  for  the  nonhydrostatic 
case.  The  numerical  models  have,  however,  been  devised  in  such 
a  way  that  more  eleiborate  dilatation  models  —  accounting,  for 
example,  for  the  influence  of  the  mean  pressure  —  can  be  imple¬ 
mented  in  a  straightforward  manner. 

Model  test  experiments  were  then  reviewed  in  an 
attempt  to  validate  the  improved  dilatation  model.  The  review 
proved,  however,  to  be  inconclusive  because: 

1.  Many  tests  have  been  performed  using  a  very  low 
friction  angle  (2.5  deg)  rock  simulant,  which, 
because  it  hardly  dilates,  can  never  provide  a 
clear  differentiation  between  full  dilatation  and 
zero  dilatation  models. 

2 .  Nvimerical  simulation  of  tests  conducted  with  a 
high- friction  rock  simulant  demonstrates  that  the 
observed  closure  is  generally  underpredicted  with 
a  linear  Mohr-Coulomb  material  even  if  a  full- 
dilatation  model  is  assvuned.  This  implies  that  a 
simple  linear  Mohr-Coulomb  criterion  is  not 
sufficient  to  describe  the  behavior  of  rock 
during  failure. 

Although  the  elastoplastic  models  developed  in  the 
course  of  this  investigation  are  based  on  a  relatively  simple 
constitutive  law,  they  nonetheless  represent  a  significant 
improvement  over  previous  analytical  models.  These  models  are, 
however,  best  used  for  parametric  analyses  and/or  to  delimit  the 
conditions  for  which  a  more  sophisticated  (but  costly)  finite 
element  analysis  is  warranted.  In  that  regard,  design  charts 
similar  to  those  developed  during  a  previous  investigation 
should  be  devised  (AA,  1983).  Such  charts  would  enhance  the 
practical  usefulness  of  the  models  developed  in  this  study. 
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APPENDIX  A 


ELASTOPLASTIC  MODEL  OF  A  DEEP  TUNNEL  FOR  A 
ROCK  WITH  VARIABLE  DILATANCY 

A-1  INTRODUCTION. 

The  hydrostatic  model  of  a  deep  tunnel  has  been  the 
subject  of  so  many  papers  (see  Brown  et  al.,  1983,  for  an 
exhaustive  list  of  references),  that  it  might  appear  unnecessary 
to  devote  yet  another  to  the  subject.  A  review  of  the  existing 
models  reveals,  however,  that  in  most  accounts  the  dilatancy  of 
the  rock  —  defined  as  the  rate  of  increase  of  the  inelastic 
volume  change  with  the  plastic  shear  strain  —  is  assumed  con¬ 
stant.  This  assumption  may  be  responsible  for  unrealistic 
prediction  of  tunnel  closure,  since  it  does  not  set  any  bound  on 
the  volume  increase  that  the  material  can  experience.  In  some 
investigations,  a  variable  dilatancy  has  been  implemented,  but 
calculation  of  the  tunnel  closure  is  then  based  on  an  approxi¬ 
mate  solution  method. 

The  objective  of  this  paper  is  to  present  a  rigorous 
solution  of  the  tunnel  closure,  for  a  general  class  of  materials 
characterized  by  a  Mohr-Coulomb  yield  envelope  and  a  plastic 
dilatation,  which  may  be  an  arbitrary  function  of  the  stress  and 
the  plastic  shear  distortion. 

A-2  THE  HYDROSTATIC  MODEL. 

Consider  the  plane  strain  model  of  a  cylindrical 

tunnel  of  radius  a,  driven  in  a  homogene  us  and  isotropic  rock 

mass  (see  figure  27).  A  far-field  stress  of  magnitude  P°  acts 

at  infinity  (it  is  assumed  that  the  gravity  force  can  be 

ignored).  Excavation  unloading  of  the  prestressed  rock  mass  is 

simulated  by  a  monotonic  decrease  of  the  internal  pressure  p, 

from  an  initial  value  P  to  zero.  The  rock  is  assumed  to  behave 

o 

as  an  elastoplastic  material  with  a  linear  Mohr-Coulomb 
envelope: 

^3  =  ^p  ^1  -  <3  (30) 
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in  which  q  is  the  unconfined  compressive  strength,  Kp  the  pas¬ 
sive  coefficient,  a  function  of  the  friction  angle,  and  tension 
and  extension  are  taken  as  positive.  We  restrict  consideration 
to  cases  where  the  out-of-plane  principal  stress  is  the  strict 
intermediate  stress  X2  in  regions  of  plastic  deformation,  so 
that  the  normal  elastic  strain,  perpendicular  to  the  plane  of 
deformation,  vanishes  everywhere.  (It  can  be  proven  that  satis¬ 
faction  of  the  inequality  (Kp  +  1)  >  l/\>  represents  a  sufficient 
condition  to  that  effect  (Detournay,  1985).) 

We  seek  to  determine  the  stress  and  displacement 
fields  in  terms  of  the  radial  coordinate  r  (the  problem  is 
axisymmetric)  and  the  internal  pressure  p,  for  a  general  class 
of  materials  characterized  by  a  plastic  dilatation  function  of 
the  stress  and  the  acciimulated  plastic  shear  strain. 

In  the  early  stage  of  unloading,  the  rock  around  the 
tunnel  remains  elastic;  however,  provided  that 
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with  P°  >  f  (31) 


the  problem  is  characterized  by  the  existence  of  a  plastic  zone 
(a  <  r  <  aR^)  surrounded  by  an  infinite  elastic  region 


surrounded 


(aR^  <  r).  Since  the  problem  is  statically  determinate,  the 
normalized  radius  R^  of  the  elastoplastic  interface  is  the  only 
function  of  the  stress  boundary  conditions  p  and  P°  and  the 
yield  parameters  q  and  Kp  (e.g.,  Salencon,  1969): 


^o  = 


P°  + 


(K„  -  1) 


K^  +  1  p  + 


(Kp  -  1) 


l/(Kp-l) 


Since  R^  is  a  monotonic  function  of  p,  it  can  be  used  as  a 
kinematic  parameter  instead  of  p  (at  least  beyond  the  elastic 
limit);  in  other  words  any  mechanical  field  (such  as  stress, 
strain,  displacement)  can  be  viewed  as  a  dual  function  of  r  and 
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The  stress  field  in  the  plane  is  independent  of  the 
flow  rule  of  the  material.  In  the  cylindrical  coordinates 
system  (r,e)  the  stresses  are  given  by  (e.g.,  Salencon,  1969): 


Plastic 


Elastic 


+  - -  V  +1 

^  Kp  -  1  ^p  ^  ^ 


+  *  V  +1 

^  Kp  -  1  ^p  ^ 


^re  =  0 


V* 


a  g  r  s  aR, 


'r  =  -  f'’  -  S? 


T  =  -  P 


•  ■:  (=t) 
■  ■■  (i;) 


r  s  aR. 


I^Xre  =  0  (34) 

in  which  the  symbol  S°  denotes  the  limiting  value  of  the  stress 
deviatoric  at  infinity: 


s°  - 

"  K  +  1 
P 


/po  +  _g — \ 

+  M  -  1/ 


within  the  elastic  region  the  induced  displacement  field  is 
given  by: 


u  =  -  ^  R_ 


r  ^  aR. 


A-3  DISPLACEMENT  IN  THE  PLASTIC  REGION. 

Within  the  plastic  region  (a  ^  r  ^  displace¬ 

ment  cannot  be  described  in  simple  form  except  for  the  case  of 
constant  dilatancy  angle.  Instead,  it  is  necessary  to  solve 
numerically  the  differential  equation  that  is  developed  in  the 
following  text. 

Since  is  used  as  the  kinematic  parameter,  the  rate 
of  change  of  a  mechanical  quantity  is  defined  as  its  partial 
derivative  with  respect  to  i  .  The  velocity  v  is  thus  defined 
by: 

V  =  ,  (37) 

o 

and  the  strain  rate  by: 


3R. 


~  aR. 


The  strain  rates  are  also  related  to  the  velocity  by: 
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In  the  plastic  region,  the  strain  (and  strain  rate) 
consists  of  an  elastic  and  a  plastic  part.  The  elastic  part  can 
be  expressed  in  terms  of  the  stress  (and  stress  rate)  by  means 
of  Hooke’s  law.  The  plastic  part  of  the  strain  rate  tensor  is 
controlled  by  the  flow  rule 

iP  =  -  K*  iP  (40) 

where  the  tangent  dilatancy  factor  K*  can  be  an  arbitrary  func¬ 
tion  of  the  stress  (generally  the  mean  pressure)  and  the  accumu¬ 


lated  plastic  shear  strain,  y  =  ”  £0-  The  flow  rule  (equa¬ 
tion  40)  can  then  be  rewritten  in  terms  of  the  total  and  elastic 
strain  rate: 


=  rP  - 


•  •  •  A  •A 

^r 


The  elastic  strain  rates  e. 


in  the  plastic  zone  can  be 


computed  from  Hooke's  law  and  the  partial  derivative  of  and 
Tq,  given  by  equation  (33),  with  respect  to  R^.  Specifically: 


S?  k. 


2G  R  (aR  ) 


V  P 


in  which 


=  (K  -  1)  (K*  -  1)  +  (1  -  2v)  (K  +  1)  (K*  +  1) 


The  differential  equation  for  the  velocity  field  in  the  plastic 
zone  can  then  be  deduced  from  equation  (41)  using  equations  (39) 
and  ( 42 ) , 


-o  .  K  -1 

^  +  K*  ^  ^ 

ar  ^r-2GR^(aRj 


a  g  r  ^  aR, 


The  velocity  field  in  the  plastic  zone  can  thus  be  calculated 
from  equation  (44)  using  as  boundary  condition,  the  value  of 
velocity  on  the  elastoplastic  interface.  This  boundary  value, 
which  is  obtained  by  differentiating  equation  (36)  with  respect 
to  R^,  and  setting  r  =  aR^  is  equal  to 

< 

V  =  -  2  ^  (45) 


Once  the  velocity  field  has  been  calculated  from  equations  (44) 
and  (45),  the  displacement  is  obtained  by  integrating  v(r,RQ) 
with  respect  to  R^.  However,  calculation  of  the  displacement 
can  be  dramatically  simplified  by  noting  that  the  displacement 
field  (like  the  stress  and  velocity)  is  actually  only  a  function 
of  the  dimensionless  ratio  r/aR^.  For  the  purpose  of  demon¬ 
stration  it  is  convenient  to  introduce  the  unit  plane  (p,9) 
defined  by  the  affine  transformation 


In  the  unit  plane,  the  circle  p  =  1  separates  an  internal  plas¬ 
tic  region  from  an  external  elastic  region  (see  figure  28). 
Using  the  transformation  (equation  46),  the  differential  equa¬ 
tion  for  the  velocity  field  becomes: 


^  +  K*  3^  -  X 

8p  p  p  * 


P  ^  1 


(47) 


in  which  v  stands  for  the  velocity  normalized  by  the  character¬ 
istic  length  L  defined  as: 

o 


as' 


L  = 


2G 


(48) 


Equation  (47),  subject  to  the  boundary  condition  v  =  -2  at 
p  =  1,  demonstrates  that  v  is  indeed  a  function  of  only  the 
cylindrical  coordinate  p  of  the  unit  plane.  Hence,  the  general 
form  of  the  displacement  field  is  necessarily  given  by: 


u(r,R^) 


(49) 


The  differential  equation  (47)  can  now  be  expressed  in  terms  of 
the  normalized  displacement  u,  using  the  fact  that 

v(p)  =  u(p)  -  p  u'(p)  (50) 


thus 


2~ 

p^u' 


+  K*  pu' 
P 


-  KJ  u  = 


-A.p 


(51) 


This  differential  equation  is  subject  to  the  boundary  conditions 


u(l)  =  -  1  ;  u'(l)  =  1 


(52) 


which  are  deduced  from  the  elastic  solution  (equation  36)  for 
the  displacement. 

If  the  dilatancy  factor  K*  is  assumed  constant 
(1  <  K*  <  Kp),  equation  (51)  is  an  Euler  equation  which  can  be 
solved  in  closed  form  to  yield 


u(p)  =  - 


\*  +  2K  +  2K* 

/  -(^5*1)  \ 

1  +  - — £ - 

(K*  +  1)(K  +  K*) 

P 

P 


This  equation  can  be  solved  readily  using  a  numerical  method, 
such  as  the  fourth-order  Runge-Kutta  technique  (Henrici,  1962), 
which  is  summarized  in  Attachment  A-1.  Applicable  solvers  are 
often  contained  within  the  math  library-  of  a  scientific  pocket 
calculator,  thus  making  the  solution  of  equation  (58)  straight¬ 
forward  even  with  limited  computational  resources. 

A-4  APPLICATION. 

As  a  simple  application  of  the  theory  developed  above, 
we  investigate  a  material  characterized  by  a  tangent  dilatancy 
factor  K*  that  decays  from  an  initial  value  according  to  an 
exponential  function  of  the  plastic  shear  strain  y: 

^  <-  f->  <59) 

v  p  y  * 

The  parameter  y^  can  most  usefully  be  related  to  the 
maximum  inelastic  volume  increase  A*,  by  integrating  the 
relation 


dy 


K*  -  1 

...P  ■ 

K*  +  1 


to  yield 


(60) 


K  +  1 

A*  =  V*  An  -^2 -  <^^) 

Equations  (59)  -  (61)  indicate  that  the  normalized 

displacement  field  u(p)  only  depends  on  three  dimensionless 

parameters,  K  ,  v,  and  A*,  which  is  defined  as 
P 


(62) 


Experimental  evidence  indicates  that  the  maximum  inelastic 
volume  increase  is  less  than  5  percent;  thus  A^  should  lie  in 
the  range  0-100. 


A- 5  CONCLUSIONS. 

The  preparation  of  this  paper  was  prompted  by  the  need 
to  improve  predictions  of  tunnel  closure.  The  assumption  of  a 
constant  dilatancy  angle  is  believed  to  be  unrealistic  because 
the  dilatation  should  be  a  function  of  the  plastic  strain  (dam¬ 
age)  and  the  confining  stress.  Accordingly,  a  variable  dilata¬ 
tion  model  was  sought.  Here  we  have  shown  that  (1)  the 
differential  equation  for  the  tunnel  closure  can  be  derived  in  a 
rigorous  manner  for  a  plastic  dilatation  which  is  an  arbitrary 
function  of  the  stress  and  the  plastic  shear  strain  and  (2)  that 
by  using  the  unit-plane  transformation,  the  differential  equa¬ 
tion  can  be  cast  in  a  form  which  is  well  suited  for  numerical 
resolution. 

*  *  *  *  * 


Attachment  1:  Solution  of  a  second-order  differential 
equation  by  the  fourth-order  Runge-Kutta  Method. 

Consider  a  second-order  differential  equation  of  the 
form  y"  =  f  (x,  y,  y'),  with  initial  values  of  x^,  y^,  y^.  The 
fourth-order  Runge-Kutta  method  leads  to  a  recursive  algorithm 
for  calculating  the  values  of  ^i+l  ^i+1  ~ 

the  known  values  of  y^,  y|  at  x^: 

n.i  =  n  "  I  "  ''4> 

^i+l  =  yi  *  »  (yj  t  i  (kj  t  kj  +  Hj)] 
where  the  coefficients  k^,  k^,  k^,  and  k^  are  given  by 

k^  =  hf(x^,  y^,  y|) 

k-  =  hf  /x.  +  ^  .  V.  +  ^  v'  +  ^  k-  .  v:  + 


APPENDIX  B 


DISPLACEMENT  FIELD  IN  THE  PLASTIC  ZONE 
CONSTANT  DILATANCY  ANGLE 


B-1  INTRODUCTION. 

This  analysis  of  the  closure  of  a  tunnel  in  a  cohe¬ 
sive,  frictional,  and  dilatant  medium,  under  nonhydrostatic 
loading,  is  based  on  the  elastoplastic  solution  derived  by 
Detournay  (1983)  (see  also  AA,  1983).  In  that  model,  the  exca¬ 
vation  of  the  tunnel  is  simulated  by  quasi-static  unloading  of  a 
hole  located  in  an  infinite  prestressed  plane.  The  particular 
load  path  selected  for  that  analysis  was  characterized  by  the 
fact  that,  beyond  the  elastic  limit  of  the  system,  unloading  of 
the  hole  corresponds  to  a  decrease  of  an  internal  pressure  p. 
Two  successive  stages  could  then  be  differentiated  in  the 
plastic  response  of  the  rock  system:  first,  the  development  of 
two  isolated  plastic  zones  on  either  side  of  the  hole,  and, 
then,  the  formation  of  a  unique  yield  region  around  the  hole. 
For  cases  where  the  hole  is  completely  surrounded  by  a  plastic 
region  and  for  cases  where  the  problem  is  statically  determi¬ 
nate,  the  equation  of  the  interface  is  given  in  complex  formu¬ 
lation  by: 

X  +  iy  =  aR^  ui(a)  (63) 

where 


u)(a)  =  \o 


(■•?) 


2/(Kp+l) 


=  [f  (-  a;  - 


(F  is  the  Gaussian 
hypergeometric  series ) 


K_  -  1 


6 


^p  ^ 


(66) 


2 

K_  +  1 


P  +  K  -  1 


K_  -  1 


l/(Kp-l) 


m  =  the  obliquity  of  the  stress  at  infinity 
a  =  the  tunnel  radius 

The  equation  (63)  for  interface  was  shown  to  be  asymptotically 
correct  for  small  departures  from  hydrostatic  loading;  nonethe¬ 
less,  it  provides  a  good  approximation  of  the  interface,  for 
cases  where  the  solution  is  statically  determinate. 

In  the  original  analysis,  the  tunnel  closure  was 
calculated  by  integrating  the  variation  of  the  incremental 
displacement  6u  at  the  tunnel  boundary  with  the  loading  para¬ 
meter;  the  incremental  displacement  6u  at  the  boundary  being 
calculated  by  solving  a  system  of  hyperbolic  partial  differ¬ 
ential  equations  governing  6u  in  the  plastic  zone,  using  as  a 
boundary  condition  the  value  of  6u  on  the  elastoplastic  inter¬ 
face.  The  implementation  of  a  variable  dilatancy  angle  neces¬ 
sitates,  however,  that  the  governing  differential  equations  in 
the  plastic  zone  be  expressed  in  terms  of  displacement  instead 
of  incremental  displacement.  As  a  first  step  toward  implement¬ 
ing  the  complete  mathematical  model  with  variable  dilatancy  (of 
Appendix  C),  we  describe  in  this  appendix,  the  new  model  for  the 
case  of  a  constant  dilatancy  angle.  After  giving  in  Sec¬ 
tion  B-2,  the  explicit  expression  of  the  elastic  displacement 
along  the  elastoplastic  interface,  we  detail  in  Section  B-3  the 
derivation  of  the  equations  governing  the  displacement  field  u 
in  the  plastic  zone  and  the  numerical  calculation  of  u  by  the 
method  of  characteristics. 

B-2  ELASTIC  DISPLACEMENT  AT  THE  ELASTOPLASTIC  INTERFACE. 

The  induced  stresses  in  the  infinite  elastic  region 
bounded  by  the  elastoplastic  interface  can  be  expressed  in  terms 
of  the  complex  potentials  and  of  Muskhelishvili 
(1962)  and  the  analytic  function  u)(0  which  maps  the  region 
exterior  to  the  unit  circle  in  the  parametric  plane  t  onto  the 


elastic  region  in  the  unit  plane  z'  =  z/aR^  (z  is  the  complex 
variable  x  +  iy  defined  in  the  physical  plane): 


'ix  ^  'w  =  2  S°  0,(0  *  4.^(0 


-  'L  *  zly  =  2  S° 


(68) 


*((0  +  ^.(Ol  (69) 

>'(0  ^  ^  J 


Where  S°  is  the  limiting  value  of  the  stress  deviatoric  S°  at 
infinity,  and  a  function  of  the  mean  pressure  P°: 

.1 


£  K  +  1 


P°  + 


■V  - 


(70) 


The  induced  displacement  in  the  elastic  region  can 
also  be  expressed  in  terms  of  the  complex  potentials  <t>,(C)  and 

X 

and  .  For  this  particular  elastoplastic  problem,  it 

can  be  shown  that  on  the  interface,  the  Cartesian  components  u  , 

x. 

Uy  of  the  elastic  displacement  are  given  by  (Detournay,  1983): 


SR,.  S, 


(71) 


where 


U  =  4(l-v)  ut(a)  <t>^(a)  -  (3-4v)  x(n) 


K_-l 


-  ^  ^  r*  P  (o)  -  i^(o) 
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a  =  e 


•'OQ 


u.(U  di 


(72) 


.  r 


u.'(U  di 


(73) 


71 


(74) 


K-1  r  _ nCK 

r*P  (o)  =  [u)(o)  P 


(K  -l)/2 


K_  -  1 


For  the  asymptotic  solution  of  the  elastoplastic  interface,  the 
Laurent  series  expression  of  the  analytic  functions  x(0/  '*'2(0/ 


<t>j^(C)  are: 


,  /f-\  -  .a  -  nii  F(-6;-6+j;  j+l;m2 

1^^^  ;2j  '  ^2j  “6  Ij/  _2, 


j=l  ^ 


F(-6;-6;l;m"-) 


with 


rrr  r 


^2j.l  =  -  2?^  r  E  ("k1  (i-k) 


F  (-6;-6+j-k;  j-k+l;m‘‘) 


i  (t)  =  T] 

^  ,.21-1  '  ^21-1  2i-l 


+  inAm_.  (77) 


with 


„„  _  ,K 

92j  -  ' 


F(-6 ;-6+j+k-l; j+k;m^ ) 


h,.  =  x''  .i- 


K_-3 


Kp*l 


(S4] 


u  ,  K  2 
h2  =  X  m 


\  1  /  \  K  +1 


Note  that  for  an  incompressible  frictionless  material, 
the  displacement  at  the  interface  is  given  by 


aR.  S 


where 


1 


b  2m  .  m\ 

-  (m^  +  l)  ct] 


Rq  =  exp 


h°.  -  P  .  i) 

\  2c  2/ 


B-3 

B-3.1 


=  c  (c  is  the  cohesion  of  the  material) 


CALCULATION  OF  THE  DISPLACEMENT  FIELD  IN  THE 
PLASTIC  ZONE. 

Governing  Partial  Differential  Equations  for 
the  Displacement. 


B-3. 1.1  Integration  of  the  Flow  Rule.  The  monotonic  load  path 
responsible  for  the  propagation  of  the  plastic  zone  around  the 
tunnel  ensures  that  there  is  no  rotation  of  the  rincipal  s-ress 
directions  in  the  plastic  zone.  Once  the  stresses  at  one  point 
reach  the  yield  surface  (i.e.,  the  point  becomes  plastic),  from 
then  on,  the  principal  stress  directions  remain  locked  along 


>  V  '  •  \ 

J'-jS'S' 


(81) 

"• 

. :  *  •  1  •» 

(82) 

.***«‘'**.' *'  • 

(83) 

the  radial  and  tangential  directions.  It  thus  follows  that  the 
principal  directions  of  the  incremental  plastic  strain  tensor 
are  radial  and  tangential  everywhere  in  the  plastic  zone  and 
at  any  time  during  the  mono  tonic  loading.  This  feature  of  the 
problem  allows  us  to  integrate  the  incremental  flow  rule 


6e 

6e 


P 

r 

P 

«1> 


(64) 


To  obtain  the  following  relation  between  the  plastic  strain 
increments  in  the  radial  and  tangential  directions: 


(85) 


The  integrated  flow  rule  (equation  85)  is  actually  equivalent  to 
the  following  two  equations  which  are  expressed  in  terms  of  the 
Cartesian  components  of  the  plastic  strain  tensor  ^ 


(®x  ■*'  ^y)  ~  (*^x  ”  ^y)  0 

2ePy  cos  2(t>  -  sin  2(|)  =  0 


(86) 

(87) 


where 


sin  <l>* 


i 

K*  +  1 
P 


(88) 


In  order  to  relate  to  the  displacement,  equations  (86)  and  (87) 
are  expressed  in  terms  of  the  total  strain,  using  the  decomposi¬ 
tion  of  the  strain  into  a  plastic  and  an  elastic  part: 
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e  P 
e  + 


(89) 


(^x  2$  -  (Ejj  -  Ey)  sin  1>* 

"  (^x  ^y)  "  fx  ■  ""y) 

( qo  ^ 


Hence : 


f  ,  w  ,  1 ;.  vy.*  V.*  njT-w^’-m^' 


«•  -■•  ’’j* 


2e^y  cos  2(t»  -  (e^  -  Cy)  sin  2$ 


=  2£j^y  cos  2(|) 


-( 


e®  -  ejj  sin  2(|> 


(91) 


The  strain  components  e  ,  e  ,  and  e  can  be  written  as  partial 

X  y  xy 

derivatives  of  the  displacement  components  u  and  u  .  Before 

^  y 

doing  so,  however,  we  derive  an  explicit  expression  for  the 
elastic  strain  in  the  plastic  zone. 

B-3.1.2  Expressions  for  the  Elastic  Strains  in  the  Plas¬ 
tic  Zone .  The  elastic  strain  in  the  plastic  zone  can  be  derived 
explicitly  as  a  function  of  the  coordinates,  using  Hooke's  law 
and  the  closed- form  expression  for  the  plastic  stresses.  Under 
the  constraint 


^z  =  0 


(92) 


which  is  assumed  to  hold,  the  elastic  stress-strain  relations  in 
the  plane  (x,y)  can  be  written  as 


®x  “  2G  ^^x  “ 
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=  2g  [(1  -  «) 


^xy  2G  ^^xy 


(93) 


Where  A  ^  denotes  variation  of  the  stress  with  respect  to  a 
reference  state  characterized  by  the  uniform  stress  Thus 

Ax. 


=  X  +  P°  -  S® 
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AXy  =  Xy  +  P°  +  S° 
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Using  equation  (94),  the  plane  stress-strain  relations  (equa¬ 
tion  93 )  can  be  rewritten  as 


*.•  1 
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(95) 


e  .  e  l-2v  _  V  .  l-2v 

X  y  2G  '  X  y'  2G 


e®  -  e®  =  ^  (x  -  T  )  - 
X  y  2G  '  X  y '  2G 


'•xy  “  2G  "^xy  '  ' 

In  terms  of  the  cylindrical  coordinates  (p,  4* )  of  the  unit- 

plane,  the  plastic  stresses  read: 


"x  V  "  ^  Kp  -  1  P 


K  -1 

tjj  -  Ty  =  2  S°  cos  241  p  P 


«  K  -1 

Txy  =  sin  20  p  P  (96) 

Using  equation  (96),  the  expressions  (95)  for  the  elastic 
stresses  in  the  plastic  zone  transform  into 


'X  ^y  \2G/ 


2(1  -  2v) 


+  1  /  K  -l\ 
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e  e 
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S°  /  K  -1  \ 

=  ^  \2  P  ^  cos  20  -  ml 

K\  Kp-^  . 

=  1^/P  20 


Where  m  is  the  obliquity  defined  as  SvS^. 

B-3.1.3  Normalized  Partial  Differential  Equations.  The  strain 

components  e  ,  £  ,  and  £  in  equations  (90)  and  (91)  are  now 
X  y  xy 

expressed  as  partial  derivatives  of  the  displacement  components 
u  and  u  .  We  will,  however,  operate  the  differentiation  in  the 

^  y 

unit  plane  (x',y')  instead  of  the  physical  plane  (x,  y).  Thus 


^x  ”  aR^  9x'  '  ^y  aR^  9y'  '  ^xy 


_±_ 

aR  \9y'  9x'  / 


since 

X  =  aR  x' 
o 

y  =  aR^  y' 

Using  equation  (98)  and  the  closed-form  expressions  (equa¬ 
tion  97)  for  the  elastic  strain,  equations  (90)  and  (91)  become: 


/  au  au  .\ 

( ^  ^  a^)  24>  - 

(ax' 

■  3^) 

sin  <l>* 

=  aR^  ^  H  (p,0) 
°  2G 


H2(P/<t>)  =  2m  sin  2$  (100) 

with 

A*  =  (Kp  -  1)(K*  -  1)  +  (1  -  2v)(Kp  +  1)(K*  +  1)  (101) 

The  system  of  partial  differential  equations  (99)  and 
the  expression  (equation  71)  for  the  displacement  at  the  elasto- 
plastic  interface  (equation  71)  represents  the  boundary  condi¬ 
tions  for  the  differential  (equation  99)  ,  indicate  that  we  can 


define  a  normalized  displacement  (u  ,  u  ),  which  is  the  only 

X  y 

function  of  the  coordinates  of  the  unit  plane; 


u  =  aR^  ^  u  (102) 

The  normalized  displacement  field  u  in  the  plastic  zone  is 
controlled  by  the  three  material  parameters  v,  K^,  K*,  and  the 
stress  obliquity  m.  Note  that  the  consistent  "normalized" 
strain  field  e,  which  is  defined  as 


=  1  +  Vh] 

i]  2 


(103) 


is  related  to  the  physical  strain  field  e  by 


e  =  ^  e  (104) 

We  now  have  to  calculate  the  normalized  displacement  components 

u  ,  u  in  the  plastic  region  of  the  unit  plane,  by  solving  the 
X  y 

following  set  of  partial  differential  equations: 


sin  4**  =  ( p  ,  0  ) 


3U 

au  \ 

+  —1 
ax'  ay'/ 
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la^ 

du 
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X  +  -JL 
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y 

sin  20  =  H2(p,0) 


with  the  boundary  conditions 


U  =  Re  [U®]  ;  u  =  Im  [U®] 

A  y 


(105) 


(106) 


along  the  curve  r '  ,  which  is  the  image  of  the  elastoplastic 
interface  in  the  unit  plane.  As  shown  in  the  next  section,  the 
system  of  equations  (105)  is  hyperbolic;  it  can  therefore  be 
solved  by  the  method  of  characteristics. 


,  ■ .... u  . ,  .  UJ  !  ■  t!  ■  .  ■  I .  L  , ,  ^  ,  .  „  ■  1^^  ^,.  ,J,  ,. 


B-3.2 


Differential  Equations  along  the  Characteristics. 


B-3.2.1  Normal  Form.  The  partial  differential  equations 
(105),  together  with  the  expressions  for  the  differentials  du^ 
and  dUy  in  terms  of  the  partial  derivatives  i.e.. 


^^x  =  ^ 


au. 


dx 


ay' 


dy' 


9%  ^'^v 


(107) 


can  be  used  to  calculate  the  first  partial  derivatives  of  u^  and 

u  at  a  point  (x',  y' )  at  which  the  differentials  du  and  du 
y  ^  y 

are  known  in  a  given  direction  dy'/dx'.  Rewriting  equa¬ 
tions  (105)  and  (107)  as  a  system  of  four  equations  in  the 

unknown  first  partial  derivatives  of  u  and  u  we  obtain: 

X  y 


cos  2if  -  sin  ♦* 


-sin  2<|i 


dx' 


cos  2^  cos  241 


dy' 


cos  2o  +  sin  ** 


sin  241 


dx' 


dy' 


3x' 


ay 


3u 

J  J 


••  *.  •, 


If  the  system  of  equations  (108)  is  hyperbolic,  there  exist  two 
real  characteristic  directions  dy'/dx',  for  which  the  determi¬ 
nant  D  of  the  system  (equation  108)  vanishes.  (Along  the 
characteristics,  the  first  derivatives  of  u  and  u  cannot  be 

^  ^  ^  y 

determined  from  the  differential  du  and  du  . )  The  vanishina  of 

x  y 

the  determinant  D  leads  to  the  quadratic  equation 

(cos  20  -  sin  d**)  -  2  sin  20 

-  (cos  20  +  sin  d>*)  =  0 

(109) 
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H-(p,0) 

du 

X 

>  (108) 

*•  -*v 1 

du 

.  y  J 

.  *  w  • 


-j 


V  ^ 


-1 


This  quadratic  equation  has  2  real  roots: 


^  -  tan  ( (|»  +  e ) 


(110) 


where 


_  n  ^  <|)* 
^  ~  4  2 


(111) 


(The  upper  sign  refers  to  the  a-,  the  lower  to  the  p-character- 
istic,  see  figure  29).  The  system  of  equations  (105)  is  there¬ 
fore  hyperbolic. 

The  differential  equation  along  the  characteristics 
are  determined  by  specifying  that  one  of  the  determinants  D ^ , 
obtained  from  D,  by  replacing  its  jth  column  by  the  column  of 
the  right  members,  is  zero  (requirement  for  a  consistent  solu¬ 
tion).  Taking,  for  example,  (corresponding  to 


^3  = 


‘cos  20  -  sin  <t>* 
-  sin  20 
dx' 

0 


cos  20  H2(P,0) 

0  du„ 


Hj^(p,0)  cos  20  +  sin  <!>’ 
sin  20 


dx' 


du. 


0 

dy' 


or 


(112) 


Dg  =  Hj^(p,0)  dx'  (sin  20  dx  -  cos  20  dy) 


-  H2(p,0)  dx'^  (cos  20  +  sin  ♦*) 


+  du^  [-2  dx'  sin  20  cos  20 


+  dy'  cos  20  (cos  20  -  sin  <t>*)] 


+  dx'  dUy  cos  20  (cos  20  +  sin  <t>*) 


(113) 


Imposing  the  vanishing  of  in  the  characteristics  direction, 
yields  for  the  a -characteristics : 


Figure  29. 


i 

I 


Figure  30 


H3^(P,0) 


du  +  tan  ( 0  -  £  )  du  =  dx '  - = - 

^  2  cos^(0  -  e)  cos  20 


(114) 


for  the  p -characteristics : 


du^  +  tan  (0  +  £)  du,  =  dx' 


Hj^(p/0 ) 

2 

2  cos  (0  +  e)  cos  20 


.  V* 

/  ^  :  o 


(115) 

The  system  of  equations  (114)  and  (115)  represents  the  normal 

form  of  the  system  of  partial  differential  equations  (105);  it 

gives  the  directional  differential  of  u  and  u  along  the  char- 

X  y 

acteristics. 

B-3.2.2  An  Explicit  Finite  Difference  Scheme.  The  normal  form 
of  the  governing  equations  of  the  displacement  field  in  the 
plastic  zone  lends  itself  naturally  to  an  explicit  finite  dif¬ 
ference  scheme.  If  equations  (114)  and  (115)  are  rewritten  as 


^i  ^^x  ^i  ^^y  =  dx'  ;  i  =  1,  2 


(116) 


The  finite  difference  discretization  of  equation  (116)  is  then 
simply  given  by  (see  figure  30) 


'•.■'vW'-, 

- 

V  j  «*• 

^ 


'*i  <“x  -  (“y  -  =y> 


=  T.  (X'  -  X' );  i  =  1,  2 
3  i 


(117) 


~3  ~3 

The  displacement  (u^,  u^)  at  point  3  —  whose  coordinates  (x^, 
y^ )  are  computed  by  calculating  the  intersection  point  of  the 
tangent  to  the  characteristics  at  points  1  and  2  -  is  thus 
calculated  by  solving  the  linear  systems  of  equations 


;  VV%  \ 


W-*  r\  wv*" 


Rl  Tj^(x^-xp  ■*■  ^1  “x  ■*■  ^1  “y 

-  ►  =  •  ►  (118) 

^2  ®2  ^2  “x  ^2  ^y 

where  the  coefficients  R^,  S^,  are  evaluated  at  (x^,  Yj^ )  • 

Numerical  tests,  using  an  algorithm  based  on  this 
scheme  (Algorithm  392  of  CACM)  revealed  poor  accuracy  of  the 
solution  (the  test  case  was  the  hydrostatic  problem),  unless  a 
high  density  characteristic  mesh  was  used.  The  problem  of 
accuracy  was  caused  by  (1)  the  curvature  of  the  characteristics, 
which  leads  to  an  error  in  the  evaluation  of  the  intersection 
point  of  the  two  characteristics,  and  (2)  the  stiff  nature  of 
the  differential  equation.  Because  of  this  problem  of  accuracy, 
another  scheme  was  implemented,  that  was  based  on  expressing  the 
components  of  the  displacements  in  the  curvilinear  coordinate 
system  of  the  characteristics. 

B-3.2.3  The  Characteristics  Coordinates  (a,  P).  The  displace¬ 
ment  characteristics  are  logarithmetic  spirals,  having  the 
origin  of  the  plane  as  an  asymptotic  point.  The  equation  of  the 
two  characteristics  intersecting  at  the  point  (p^,  are  given 


P  =  Po  e 


±(<t'Q-(t>)  tan  (7t/4-<l>V2) 


(119) 


(upper  sign  a -characteristics ,  lower  sign  ^-characteristics ) . 

It  follows  from  equation  (119),  that  the  a-  and 
^-characteristics  can  be  identified  as  the  coordinate-lines  of  a 
curvilinear  coodinates  system  (a,8).  (Constant  p  and  a  coordi¬ 
nates,  respectively)  defined  as: 


«  =  -  ^K*£np+<ti 

P  =  -  ^  £n  p  -  (t>  (120) 

The  characteristic  coordinates  (a,  p)  have  been  defined  in  such 
a  way  that  the  base  vectors  e^  and  e^  are  pointing  towards  the 
asymptotic  point  (see  figure  29).  The  contravariant  components 
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and  Up  in  the  characteristic  coordinate  system  (a,  p)  are 
given  by  (see  figure  31). 


■  -*• 

1 


u^  =  -  cos  (0  -  e)  Ujj  -  sin  (0  -  e)  Uy 


u„  =  -  cos  (0  +  e)  u„  -  sin  (0  +  e)  u 
P  X  y 


(121) 


B-3.2.4  Ordinary  Differential  Equations  Along  the  Character¬ 
istics  .  The  governing  differential  equations  of  the  displacement 
field  in  the  plastic  zone  can  now  be  rewritten  in  terms  of  the 
contravariant  components  of  the  displacement.  Inverting  equa¬ 
tion  (121), 


K 


r 


cos  0*  Uj^  =  - 

sin  (0  +  e) 

sin  (0 

- 

£  ) 

u 

P 

cos  0*  Uy  =  - 

cos  (0  +  e ) 

- 

cos  (0 

- 

e  ) 

u 

P 

( 

differentiating  the 

above  equations, 

we  obtain 

cos  0*  du^  =  - 

sin  (0  +  e ) 

du 

a 

- 

cos 

(0 

+ 

e  ) 

d0 

+ 

sin  (0  -  e ) 

dUp 

+ 

cos 

(0 

- 

e  ) 

^P 

d0 

cos  0*  dUy  = 

cos  (0  +  e ) 

du 

a 

- 

sin 

(0 

+ 

e  ) 

d0 

+ 

cos(0  -  e  ) 

dSp 

+ 

sin 

(0 

- 

e  ) 

”p 

d0 

(122) 


(123) 

Substituting  and  du^  as  given  by  equation  (123)  in  the 
differential  equations  (114)  and  (115),  and  expressing  dx'  as  a 
function  of  d0  leads  to: 


du  -  (u  tan  <l>*  +  u^  sec  il>*)  d0 
a  '  a  P  I  -r 


-£di_ 


2  sin  e  cos  20 


jH^(p,0)  +  2  sin  (0  -  e)  cos  (0  -  e)  (p,0)] 


(124) 
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du_  +  (u-  tan  +  u  sec  ♦*)  d*  =  -  ■= — .  - — 

P  P  a  '  ^  2  sin  e  cos  2^ 

[Hj^(p,0)  +  2  sin  ($  +  e)  cos  +  e)  (p.<t>)] 

(125) 


Since 


dx'  =  -  r  pd<l>  along  a -characteristics 

sxn  c 

dx'  =  t  along  p-characteristics 

S  ^Xl  & 

and  noting  from  equation  (120)  that 
dof 

d(^  =  ^  along  p  =  constant 


d(|)  =  -  ^  along  a  =  constant 


the  differential  equations  (124)  and  (125)  can  finally  be 
rewritten  as 


2  hT-  -  (u  tan  +  u.  sec  ♦*)  =  t 


(126) 


8  ~  ~ 

2  -  (Up  tan  4>*  +  u^  sec  <»>*)=  tp 


(127) 


where,  after  simplification,  t^  and  tp  are  given  by 


t  =  sin  £  p  -  - j —  p 

(Kp-l)(Kp+l) 

K_  +  l 


+  (l-2v)  -  m  cos  2  (()»-£) 


(128) 


t  =  sin  £  p - - —  p  P  +  (l-2v) 

P  (K  -1)(K  +1)  *^p‘-^ 


-  m  cos  2  (0-*-£  ) 


The  cylindrical  coordinates  (p,i^)  in  the  above  equations  can 
simply  be  expressed  in  terms  of  a  and  p  by 


(130) 


B-3.3  Numerical  Calculation  of  the  Displacement 

Field  in  the  Plastic  Zone. 

The  numerical  determination  of  the  disp^lacement  field 
in  the  plastic  zone  is  carried  out  by  the  method  of  characteris¬ 
tics  (Masseau,  1899),  which  is  based  on  the  discretization  of 
the  differential  equations  (126)  and  (127). 

Consider  the  point  P,  inside  the  plastic  zone  (see 
figure  30).  The  displacement  at  P  is  controlled  by  the  values 
of  the  elastic  displacement  along  the  arc  AB  of  r,  which  is 
intercepted  by  the  two  characteristics  intersecting  at  P.  (In 
other  words,  the  domain  of  determinacy  of  the  arc  AB  is  the 
curvilinear  triangle  ABP  bourUed  by  the  two  characteristics  a 
and  p.)  In  the  method  of  characteristics,  the  displacement  at  P 
is  approximately  solved  by  first  defining  N  nodes  along  the 
noncharacteristic  arc  AB,  then  progressively  computing  the 
displacement  at  all  the  nodes  of  the  characteristic  mesh 
( located  at  the  intersection  of  the  characteristics  emerging 
from  the  initial  nodes),  using  a  discretized  form  of  the  dif¬ 
ferential  equations  (126)  and  (127). 

In  the  following,  we  derive  the  equations  needed  to 
calculate  the  displacement  at  any  node  (a^,  p^),  assuming  the 
displacements  known  at  the  two  "parent"  nodes  (a^,  p^^)  and  (a^, 
P2).  (03  =  ~  ^1’  ^ 

A  class  of  numerical  algorithms  to  calculate  the 
displacement  at  the  node  (02,  Pj^)  is  based  on  the  so-called 
4-method.  The  4"n*®thod  relies  on  two  assumptions: 

1.  The  contravariant  components  (u^,  u^ )  vary 
linearly  with  the  characteristic  oordinate 
between  two  adjacent  nodes  located  on  the  same 
characteristics . 
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The  differential  equations  (126)  or  (127)  hold  at 
a  certain  point  •'4*'  of  the  characteristic  arc 
defined  by  two  adjacent  nodes. 

For  example,  consider  the  point  ”4”  of  coordinates  (a,  p^^) 
located  on  the  characteristic  arc  defined  by  the  two  end  nodes 
(Oj^/  Pj^)  and  (02,  Pj^)-  The  value  of  i,  which  must  be  in  the 
range  (0,1)  is  given  by 


i  = 


0-0. 


“2  -  “l 


(132) 


The  displacement  at  point  4  is  given  by 


Ci 

Q 

II 

(1 

-  4)  Si 

^ '  a 

+  4 

0 

II 

oa 

(1 

-  s»  “J 

+  4 

G3 

0 

(133) 

On 

the 

basis 

of 

equation  ( 133 ) 

,  the 

differential 

126, 

which  is  assumed  to  hold  at 

point 

"4",  becomes: 

^  K  ‘  -  <"2.‘  “i>  **  ^o  ^ 


+  sec  <t>* 


[|  .  (i-t)  3j] 


+  t 


=  0 


(134) 


where  t  is  used  to  denote  the  value  of  t  at  point  "4"-  The 
discretization  of  equation  127  can  be  carried  out  in  a  simi¬ 
lar  way: 

2  [“s  -  -  ^2'  1'^“  ♦*  [«  “I  *  Up] 

+  sec  4.*  [4  +  (1-4)  5^]  +  tpj  =  0 

(135) 

The  above  two  equations  can  be  written  in  matricial  form  as: 


(136) 


1^00 

^op 

[5^1 

0 

•  =  - 

•  ^ 
p 

a 

_°pa 

°PP, 

L  pJ 

KJ 
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where 


2/Aa  -  4  tan  ♦* 

-  4  sec 

-  4  sec  ♦* 

2/Ap  -  4  tan  <l>* 

+  tan  <!>*)  5^  +  (D^p  +  sec  4>*)  uj  +  t^ 

(Dpp  +  tan  «t.*)  S2  +  (Dp^  +  sec  ♦*)  +  tp 

with 

Aa  =  a2  - 
AP  =  Pi  -  P2 

A  class  of  algorithms  can  thus  be  generated  depending  on  the 
value  of  4/  from  full  explicit  (4  =  0)  to  fully  implicit 
(4  =  !)•  Some  parametric  investigations  confirmed  that  the 
central-difference  method  (4  =  0.5)  provides  the  most  accurate 
scheme . 

~3 

For  the  central -difference  method,  the  components  (u  , 

^  3  ^ 

Up)  of  the  displacement  at  node  k  are  given  explicitly  by 

-I  =  -  5^  *  (l  -  tan  f ^  *  =ec  ♦*  ^ 

"p  =  ■  "p  *  (l  '  tan  ♦*  1^)  ^  +  aat  ♦*  |^  ^  (137) 

where 

'>«  =  (=«  t  .  tan  »*  -  52)  .  t„ 

‘>p  '  (“J  t  tp)  t  tan  **  f  (S2  -  Sj)  .  tp  |i 

D  =  (1  -  tan  ♦*  1  -  tan  0* 


Aa  At 


1  £. 


2 

(sec  <l>* 


tan^  <!>*) 


APPENDIX  C 


DISPLACEMENT  FIELD  IN  THE  PLASTIC  ZONE-VARIABLE  DILATANCY 


C- 1  INTRODUCT I ON . 

In  this  appendix,  we  develop  the  theoretical  basis  of 
a  numerical  algorithm  to  calculate  closure  of  a  circular  tunnel 
subject  to  a  nonhydrostatic  far-field  stress,  for  a  class  of 
materials  characterized  by  a  variadsle  dilatancy.  The  following 
analysis  closely  parallels  the  one  outlined  in  Appendix  B,  for  a 
material  with  a  constant  dilatancy. 

C-2  VARIABLE  DILATANCY  FACTOR. 

The  dilatancy  factor  K*,  which  is  defined  as 


is  assumed  to  decay  exponentially  with  the  accumulated  plastic 
shear  strain  y,  from  an  initial  value  Rpt 

KJ  =  1  +  (Kp  -  1)  (139) 


The  flow  rule  (equation  139)  is  associated  at  the 
elastic  limit  (K*  =  Kp  if  y  =  0);  but  as  the  material  is  yield¬ 
ing,  the  rate  of  increase  of  the  inelastic  dilatation  A  with  y 
progressively  decreases  so  as  to  eventually  vanish,  when  the 
■maximum  dilatation  A*  is  reached.  The  parameter  y*  in  equa¬ 
tion  (139)  can  be  related  to  by  integrating 


dy  KJ  +  1 


(140) 


to  yield 


90 


where 


with 


%  (p.0)  =  - 


(<S,  - 1)  *  1) 


le-i 

p  ^  cos  2^ 


K*  -  1 


+  2(l-2u)  ^ - r  cos  2^  +  2m  - 

\  ~  ^  Rj  +  1 


^2  ~  2in  sin  2^ 


=  <Kp  -  1)(K*  -  1)  +  (1  -  2v)(Kp  +  1)(K*  +  1) 


(146) 


(147) 


y^:yy.y 

-w*  -  •  • 

*>  •  • 


■^v 


Despite  their  similarity,  an  important  difference  exists  between 
the  two  sets  of  equations  (99)  and  (145).  In  the  case  of  a 
constant  dilatancy,  the  normalized  displacement  field  u  in  the 
plastic  domain  of  the  unit-plane  depends  on  the  material  para¬ 
meters  Kp,  K*,  y,  and  the  stress  obliquity  m;  hence  the  same 
normalized  displacement  field  holds  for  any  shear  elastic  modu¬ 
lus  and/or  any  stress  at  infinity  characterized  by  the  same 
obliquity  m.  In  the  case  of  a  variable  dilatancy  angle,  how¬ 
ever,  the  governing  equations  for  u  depend  on  the  ratio  y/^*  (by 
virtue  of  the  law  of  variation  of  K*)  besides  the  pareuneters  K  , 

c' 

V,  m.  Since 


^  =  2G  ^ 


(148) 


where 


au  u  1  au. 

’  ap  p  p  90 


Kp-1 


(149) 


the  coefficients  of  the  system  of  ecjuation  (145)  are  actually 
functions  of  the  ratio  y/A*,  where  is  given  by 


.••/v 


Consequently,  the  normalized  plastic  displacement  field  u 
depends  now  on  the  dimensionless  parcuneters  m,  K^,  v,  and 
(The  same  conclusions  would  hold  for  other  forms  of  variation  of 
K*.)  Hence,  some  of  the  properties  of  the  solution  for  a  con¬ 
stant  dilatancy  angle  (e.g.,  independence  of  u  on  S°,  G)  do  not 
hold  anymore. 

C-3.2  Differential  Equations  along  the  Characteristics. 

Derivation  of  the  differential  equations  along  the 
characteristics  is  identical  to  the  procedure  detailed  in  Sec¬ 
tion  B-3  of  Appendix  B,  but  for  the  substitution  of  <l>*  by 
The  normal  form  of  the  differential  equations  is 

du  +  tan  ( 0  -  £  )  du  = 

X  y 

+ 


du  +  tan  ((>+£)  du  = 
X  y 

+ 


and  the  characteristic  directions 

^  =  tan  (0  +  £)  (152) 


dx 


Hj^(p»'t') 


2  cos  ((>  -  e)  cos  20 


along  the  or -characteristics 

(151) 

H^(P,<)») 


dx 


2  cos  (0  +  e)  cos  20 


tan  ( 0  +  e ) 

^  SSs  2<^ 

along  the  p-characteristics 

are  given  by 


(upper  sign  for  a-,  lower  sign  for  p-characteristics).  In 
equations  (151)  and  (152),  the  symbol  e  denotes  the  inclination 
of  the  characteristics  on  the  radial  direction: 


Since  the  inclination  e  is  a  function  of  the  solution 
(i.e.,  the  displacement  u),  it  is  not  anymore  possible  to  define 
explicitly  the  characteristic  coordinates  (a,p).  However,  the 
formal  form  of  the  differential  equations  (126)  and  (127)  still 
holds ;  i.e.. 


^ . 
da 

tan 

'^p 

sec  $*)  =  t^ 

(154a) 

dp 

tan 

+  u 

a 

sec  ♦*)  =  tp 

(154b) 

Provided  that  da  and  dp  are  defined  by 

d^  =  ft  .  (155) 

Indeed,  the  relationship  between  the  increments  dp  and  di|)  - 
characterizing  the  variation  of  the  cylindrical  coordinates 
between  2  points  infinitesimally  close  on  the  same  characteris¬ 
tics  (see  figure  32)  -  is  given  by  da  =  0  for  the  p -characteris¬ 
tic,  and  by  dp  =  0  for  the  a -characteristic.  Note  that  in 

equations  (154)  the  curvilinear  components  u^  and  Up  are  given 

in  terras  of  u  ,  u  ,  by  equation  (121),  with  e  replaced  by  e,  and 
X  y 

the  values  of  t^  and  tp  are  obtained  from  equations  ( 128 )  and 
(129)  respectively  with  e,  K*,  replaced  by  e,  K*,  \*. 

C-4  NUMERICAL  SOLUTION  OF  THE  DISPLACEMENT  FIELD. 

C-4.1  Preeunble. 

The  system  of  differential  equations  (154)  will  be 

solved  by  the  method  of  characteristics.  A  sequence  of  N  nodes 
is  selected  along  the  elastoplastic  interface,  that  defines  a 
fan  of  characteristics  in  the  plastic  domain.  The  displacement 
at  the  initial  nodes  on  r  are  calculated  from  the  solution  of 
the  displacement  field  in  the  elastic  domain,  while  the  dis¬ 
placement  at  the  nodes  of  the  characteristic  mesh  is  progres¬ 
sively  computed  by  moving  away  from  r ,  using  the  discretized 
form  of  the  differential  equations.  These  discretized  equations 
will  be  derived  for  a  class  of  methods  (the  t-method),  which 


do,  =  ^  ; 

0 


provides  a  full  spectrum  between  fully  explicit  to  fully 
implicit. 

Several  complications  arise  in  the  calculation  of  the 
plastic  displacement  field,  that  are  introduced  by  the  varicible 
character  of  the  material  dilatancy.  First,  the  position  of  the 
nodes  of  the  characteristic  mesh  cannot,  as  previously,  be 
calculated  prior  to  the  displacement  of  these  nodes.  Instead, 
they  are  an  integral  part  of  the  solution  and  must  be  calculated 
concurrently  to  the  displacement  at  the  nodes.  Second,  the 
differential  equations  are  nonlinear,  and  so  are  their  dis¬ 
cretized  form  (except  for  the  fully  explicit  case).  As  a  con¬ 
sequence,  the  calculation  of  the  displacement  at  a  new  node,  and 
the  position  of  this  node  requires  an  iterative  computational 
procedure.  The  following  sections  detail  the  basis  of  a  numeri¬ 
cal  algorithm  to  calculate  the  displacement  field  in  the  plastic 
zone. 

C-4.2  Discretized  Equations. 

Consider  two  close  points  1  and  2  of  the  plastic  zone, 

at  which  the  displacement  is  known  (see  figure  33).  Let  p^, 

denote  the  cylindrical  coordinates  of  point  i  (i  =  1,2)  and  u^, 

uj  the  cylindrical  components  of  the  known  displacement  at  those 

points.  For  the  sake  of  definiteness,  it  is  assumed  that 

1^2  ^  We  need  to  calculate  the  coordinates  (p^,  0^) 

point  3,  the  intersection  of  the  or-characteristic  through 

point  1  and  the  p -characteristic  through  point  2,  and  the  dis- 
~3  ~3 

placement  u^ )  at  that  point. 

First  introduce  the  characteristic  "coordinates"  a,  p 
which  are  only  valid  on  the  two  arc  segments  13  and  23 


on  23 : 


(*p)f 


p"  =  - 


(‘^)i 


jKn  p  + 


£n  p  -  41 


(157) 


where  (K*)„  and  (K*)„  represent  an  "average"  value  of  the  secant 
dilatancy  factor  along  segments  13  and  23,  respectively 
(actually  the  value  of  K*  at  point 

The  characteristic  "coordinates"  (a^,  pp  of  point  1 
are  given  by  equation  (156),  with  p  and  0  substituted  by  p^^  and 
0^.  Similarly,  the  "coordinates"  (a",  pp  of  point  2  are 
obtained  from  equation  157. 

The  segment  of  characteristic  13  is  characterized  by 


p'  =  p' 


(158) 


'f', 

•  V  V  ■ 


and  12  by 


a"  =  a'^ 


(159) 


As  for  the  case  of  a  constant  dilatancy  angle,  we  can  define  the 

curvilinear  components  u^  and  u^  of  the  displacement,  which  are 

related  to  the  cylindrical  components  u  ,  u^  by 

P  0 


u  =  - 
a 


^p 


■  u 
K*  +  1  P 
P 


.  -.=■■■  u  + 
K*  +  1  P 
P 


'K*  +  1 
P 


(160) 


K*  +  1 
P 


where  K*  is  either  (K*)  or  (K*)„,  depending  on  whether  the 
P  pa  P  p 

point  is  on  the  13  or  23  arc  segments. 

The  numerical  procedure  to  calculate  p^,  03,  u^,  u^ 

relies  on  the  following  assumptions: 

1.  The  curvilinear  components  u^  and  Up  vary 

linearly  with  a'  on  the  arc  13,  and  with  p"  on 

the  arc  23. 


2. 


The  differential  equation  (154a)applies  at  point 
”4''  (0  5  4^1)  of  segment  13,  with  coordinates 
{4a^  +  (1  ~  4)  «2'  equation  (154b)  at 

the  coordinates  (a",  4  P2  (i  "  O  on  the 

arc  23. 

3.  At  point  3,  =  Og  =  0(3  and  P3  =  P3  s  P3,  which 

is  only  correct  if  (K*)^  and  (K*)^  correspond  to 

pa  p  p 

the  value  of  the  secant  dilatancy  factor  at 

3  3 

point  3,  and  if  u^  (and  Up )  is  equal  in  both 
characteristic  systems  («',  P')  and  (a",  p"). 

4.  The  "average"  secant  dilatancy  factor  (K*)^  is 
equal  to  the  value  of  Kp  at  point  "4"  on  seg¬ 
ment  13;  similarly  (K*)-  is  the  value  of  K*  at 

P  p  P 

point  "4"  on  segment  23. 

It  follows  from  the  assumptions  that  the  discretized  form  of  the 
two  differential  equations  (154)  are 

{h  -  5  *s)  “o  - «  ===  “3 

-  5  sec  ♦*  .(Ij  -  {  tan  ♦*)  ^ 

=  (fj  t  (1  -  {)  tan  «*)  52.  (1  .  j,  sec  .  t^ 


where 


(161) 


t  and  t  =  t  and  t„  calculated  at  point  "4"  on  seg- 
“  ^  ments  island  23,  respectively 

$*  and  =  the  secant  dilatancy  angle  at  point  "4"  on 
“  ^  segments  13  and  23,  respectively 


the  discretized  equations  (163)  can  be  rewritten  in  matricial 
form: 


with 


^pa  °ppj  l^p 


^ap  =  ■  ^  H 

Dpa  =  -  ^  ^P 


(162) 


^pp 

=  2/AP 

-  i 

tan  ♦* 
p 

^^gg 

+  tan 

5*) 

*  '“op 

+  sec  $*) 

+  t 

=  <°pp 

+  tan 

$*) 

5p  *  ">p„ 

+  sec  ♦?) 
p 

S2 

g 

+  t 

^3  ^3 

The  components  u  and  of  the  displacement  at  point  3  are  then 
a  p 


given  by 


u3  _(  ^pp  “  ^p  ^gp) 

a  \  D  / 

53  ,(  °»a  -  °p.) 


(163) 


with 


^  "  °aa  °pp  "  °ap  °pa 

~3  ~3 

Finally,  the  cylindrical  components  (u^,  )  of  the  displacement 


at  point  3  read 

“r  =  - 1  (V<¥^  * 

l/(K*)p  +  l 


)(=^  *  'i) 


'(K*)  +1 
p  g 

V(K*) 

p'g 


(164) 
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while  the  position  of  point  3  is  given  by 
Po  =  exp 


V<^  “2  ^ 


(165) 


C-4.3  Iterative  Procedure. 

The  quantities  (K^)^  and  (K^)q,  and  the  related  con- 

p  a  p  p 

stants  which  appear  in  the  systems  of  equations  (162),  are  not 
known  beforehand  since  they  are  functions  of  the  plastic  shear 
distortion  y  at  point  "4".  An  exception,  however,  is  the  fully 
explicit  case  for  which  point  "I”  on  segments  13  and  23  cor¬ 
responds  to  the  known  points  1  and  2,  respectively.  Thus,  but 
for  the  case  4  =  0,  the  system  of  equations  162  is  nonlinear 
and  has  to  be  solved  iteratively.  The  iterative  procedure 
consists  of  taking  as  a  first  approximation  of  (Kp)^  and  (Kp)^, 
the  values  of  Kp  at  1  and  2,  solving  the  systems  of  equa¬ 
tions  162,  calculating  the  positions  and  displacements  of 

k  ^k 

points  ”4"  to  determine  new  approximations  for  (Kp)^  and  (Kp)p* 
and  iterating  until  satisfactory  convergence  is  achieved. 

C-4.4  Calculation  of  the  Plastic  Distortion  y- 

The  plastic  distortion  y  can  be  calculated  from  the 

flow  rule 

A  =  sin  y  ( 166 ) 

and  the  knowledge  of  the  plastic  normal  strain  in  a  non¬ 
characteristic  direction.  Indeed,  assume  that  is  known  in  a 
direction  which  is  inclined  by  an  angle  n  on  the  x-axis;  thus 

eP  =  1  (eP  +  eP)+  i  cos  2  (n  -  «D)  -  ej  (167) 

Using  equations  (166)  and  (167),  we  obtain  for  y 


99 


(168) 


^  •• 


sin  ♦*  +  cos2(r|-<t>) 

If  n  defines  a  characteristic  direction,  cos  2{n-(|»)  =  -  sin  '!>* 
and  \  is  undefined  from  equation  (168)  (e^  vanishes  in  the 
characteristic  direction). 

Let  us  now  calculate  the  approximate  normal  strain  e 
between  two  adjacent  points  A  and  B  at  which  the  displacement  is 
known  (see  figure  34).  The  angle  n»  which  gives  the  inclination 
of  the  segment  AB  on  the  x-axis,  and  the  distance  AL  between  A 
and  B  are  given  by 


n  =  arc  tan 


yB  -  ^A 

*B  -  *A 


AL 


=  V  <’‘b  - 


(169) 


Let  u-o  denote  the  displacement  at  point  A  in  the  n -direction 

^  A 

and  u^  the  displacement  at  B,  also  in  the  n "direction,  u^  and 

u^  are  related  to  the  cylindrical  components  of  the  displace¬ 
ment  at  A  and  B  by 


Uab  =  cos  (n-<|i)  Up  +  sin  (n-4>)  u^  (170) 

The  average  extension  e  between  A  and  B  is  thus  given  by 


e 


AL 


(171) 


The  average  strain  e ,  represents  an  approximation  of  the  normal 
strain  e  in  the  n "direction  at  both  points  A  and  B. 

The  elastic  past  e®  of  the  normal  strain  in  the  direc¬ 
tion  n  is  given  in  terms  of  the  Cartesian  components  of  the 
elastic  strain  tensor  by 


~e 

e 


cos  2a  +  e® 
xy 


sin  2a 


(172) 


Since 


5  i^x  ~  ^y)  “  P  ^  2^  -  m 

~e  . 

Cjjy  =  P  Sin  2$ 


(173) 


the  expression  for  e  becomes 


e®  =  (1  -  2v 


+  1  /  K„-l\ 

>  K^^T  (i  -  0  ) 


K  -1 

+  cos  2(n-0)P*^  -™  cos  2r) 

The  plastic  part  is  then  determined  by  the  difference  e  -  e’ 
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